Implement k-point grid

Author: vitartas

app/quoll.jl

@@ -68,12 +68,16 @@ for idir in my_idirs
 
     if Quoll.Parser.requires_kpoint_grid(params)
         @info "Initialising k-point grid"
+
+        # Assuming all operators in the directory are related to the same atoms
+        atoms = Quoll.get_atoms(first(operators))
+        kgrid = get_kgrid(atoms, params)
     end
 
     if !isnothing(params.output)
         @info "Converting operators into $(params.output.format) format"
 
-        operators = convert_operators(operators, params.output.format; hermitian = params.output.hermitian)
+        operators = convert_operators(operators, params.output.format, hermitian = params.output.hermitian)
 
         @info "Writing operators"
 

examples/SiC/input_file.toml

@@ -11,7 +11,7 @@ directory = "SiC_DeepH_converted"
 # Should also allow different types,
 # i.e. Monkhorst-Pack and even possibly different shifts
 [kpoint_grid]
-grid = [10, 10, 10]
+mesh = [10, 10, 10]
 density = 8.0
 
 [basis_projection]

src/Parser.jl

@@ -12,11 +12,12 @@ using ..Quoll
 using ..Quoll:
     AbstractOperator,
     OperatorKind,
+    allows_symmetry,
     get_avail_operatorkinds,
     get_operatorkinds,
     get_readformat,
     get_writeformat
-import ..Quoll: find_operatorkinds
+import ..Quoll: find_operatorkinds, get_kgrid
 using ..Projections
 using ..Projections: AbstractBasisProjection, get_basis_projection
 using ..Postprocessing
@@ -63,12 +64,17 @@ include("parser/symmetry.jl")
             postprocessing,
             error_metrics,
         )
+        # Check if supplied symmetry options are appropriate
+        validate_symmetry(input, basis_projection, symmetry)
+
+        # Look if there are any unresolvable clashes
         search_clashes(basis_projection, error_metrics)
+
         new(input, output, basis_projection, kpoint_grid, postprocessing, error_metrics, symmetry)
     end
 end
 
-# Methods that either require all subparameter types or QuollParams itself
+# Methods that either require subparameter types or QuollParams itself
 include("parser/methods.jl")
 
 end

src/Quoll.jl

@@ -14,15 +14,20 @@ using JSON
 using AtomsBase
 using AtomsIOPython
 using NeighbourLists
+using Spglib
 using MPI
 
-# Common tools
+# Constants
+
+include("constants.jl")
+
+# Methods
 
 export recentre
-include("common/constants.jl")
-include("common/mpitools.jl")
-include("common/atomtools.jl")
-include("common/fastlookup.jl")
+include("methods/mpitools.jl")
+include("methods/atomtools.jl")
+include("methods/fastlookup.jl")
+include("methods/symmetry.jl")
 
 # Core types and their methods
 
@@ -37,6 +42,9 @@ include("spin.jl")
 include("sparsity.jl")
 include("shconversion.jl")
 
+export get_kgrid
+include("kgrid.jl")
+
 # Operator formats and their IO routines
 # (Assuming no dependencies between different format types in each file)
 

src/common/atomtools.jl

@@ -4,11 +4,8 @@ const AbstractAnyDict{K, V} = Union{<:AbstractDict{K, V}, <:AbstractDictionary{K
 
 const SpeciesPairAnyDict{V} = AbstractAnyDict{NTuple{2, ChemicalSpecies}, V}
 const SpeciesAnyDict{V}  = AbstractAnyDict{ChemicalSpecies, V}
-# const SpeciesDict{V} = AbstractDict{ChemicalSpecies, V}
-# const SpeciesDictionary{V} = AbstractDictionary{ChemicalSpecies, V}
 
 const AtomPairAnyDict{V} = AbstractAnyDict{NTuple{2, Int}, V}
-# const AtomPairKeyedArray{T, N, AT, KT} = AbstractArray{<:KeyedArray{T, N, AT, KT}, 2}
 
 const ThreeDimSliceView{T} = SubArray{
     T, 3, Array{T, 3},

src/common/constants.jl src/constants.jlsrc/common/constants.jl renamed to src/constants.jl

@@ -0,0 +1,97 @@
+
+# TODO: move some of the functions from here to atomtools.jl
+# and some file containing symmetry operations like obtaining rotations
+# and checking if the system has inversion symmetry or not (relevant for Brillouin.jl)
+
+# TODO:
+# - Would be useful to know whether it contains symmetries or not, and which ones
+# - There needs to be a way to distribute k-points to MPI tasks
+#   (to achieve this could have an outer constructor which would take full KGrid
+#   plus my_indices and would return a smaller KGrid)
+struct KGrid{T, W}
+    kpoints::T
+    weights::W
+    time_reversal::Union{Missing, Bool}
+    crystal_symmetry::Union{Missing, Bool}
+end
+
+function get_kgrid(atoms::AbstractSystem; density = nothing, mesh = nothing,
+    shift = falses(3), time_reversal = false, crystal_symmetry = false, symprec = 1e-5)
+
+    # Get symmetry rotations
+    spglib_cell = get_spglib_cell(atoms)
+    rotations = get_symmetry_rotations(spglib_cell, crystal_symmetry = crystal_symmetry, symprec = symprec)
+
+    # Get k-point grid
+    mesh = get_recip_mesh(atoms, density, mesh)
+    kgrid_spglib = get_stabilized_reciprocal_mesh(rotations, mesh, is_shift = shift, is_time_reversal = time_reversal)
+
+    # Add the shift and wrap irreducible k-points to [-0.5, 0.5)
+    kirreds = map(Spglib.eachpoint(kgrid_spglib)) do kcoord
+        normalize_kpoint_coordinate(SVector{3}(kcoord))
+    end
+
+    # Get groups of reducible k-points
+    ir_mapping = kgrid_spglib.ir_mapping_table
+    k_all_reducible = [findall(isequal(elem), ir_mapping) for elem in unique(ir_mapping)]
+
+    # Get k-point weights
+    n_equivalent_k = length.(k_all_reducible)
+    weights = n_equivalent_k / sum(n_equivalent_k)
+
+    # Test if all reducible k-points can be reproduced with symmetry rotations
+    grid_address = kgrid_spglib.grid_address
+    symmetries_valid = check_kpoint_reduction(rotations, time_reversal, shift, mesh, grid_address, k_all_reducible, kirreds)
+
+    if !symmetries_valid
+        @warn "Reducible k-points could not be generated from the irreducible kpoints. " *
+              "This points to a bug in spglib. Defaulting to no symmetries."
+        return get_kgrid(
+            atoms, density = density, mesh = mesh, shift = shift,
+            time_reversal = false, crystal_symmetry = false, symprec = symprec
+        )
+    else
+        return KGrid(collect(eachpoint(kgrid_spglib)), weights, time_reversal, crystal_symmetry)
+    end
+end
+
+"""Bring ``k``-point coordinates into the range [-0.5, 0.5)"""
+function normalize_kpoint_coordinate(x::Real)
+    x = x - round(Int, x, RoundNearestTiesUp)
+    @assert -0.5 ≤ x < 0.5
+    x
+end
+normalize_kpoint_coordinate(k::AbstractVector) = normalize_kpoint_coordinate.(k)
+
+function check_kpoint_reduction(rotations, time_reversal, is_shift, mesh, grid_address, k_all_reducible, kirreds)
+    all_rotations = time_reversal ? rotations ∪ [SMatrix{3, 3}(-I(3))] : rotations
+    shift = is_shift .// 2
+
+    for (iks_reducible, k) in zip(k_all_reducible, kirreds), ikred in iks_reducible
+        kred = (shift .+ grid_address[ikred]) .// mesh
+        found_mapping = any(all_rotations) do W
+            # If the difference between kred and W' * k == W^{-1} * k
+            # is only integer in fractional reciprocal-space coordinates, then
+            # kred and S' * k are equivalent k-points
+            all(isinteger, kred - (W * k))
+        end
+        if !found_mapping
+            return false
+        end
+    end
+    return true
+end
+
+function get_recip_mesh(atoms::FlexibleSystem, density, mesh)
+    if isnothing(mesh)
+        real_lattice = ustrip.(stack(cell_vectors(atoms)))
+        return get_recip_mesh(real_lattice, density)
+    else
+        return mesh
+    end
+end
+
+function get_recip_mesh(real_lattice, density)
+    recip_basis_lengths = norm.(eachcol(get_recip_lattice(real_lattice)))
+    return ceil(Int, recip_basis_lengths .* density)
+end

src/kgrid.jl

@@ -0,0 +1,48 @@
+function recentre(atoms::AbstractSystem)
+    # Currently the following function does not support recentering
+    # for systems with open boundary conditions
+    @argcheck all(periodicity(atoms))
+
+    # Convert atom positions to fractional coordinates.
+    atom_frac_pos = get_frac_positions(atoms)
+
+    # Wrap fractional coordinates to [-0.5, 0.5) unit cell
+    # and shift back to [0.0, 1.0)
+    wrapped_atom_frac_pos = atom_frac_pos .- 1e-8
+    wrapped_atom_frac_pos = wrapped_atom_frac_pos .- round.(wrapped_atom_frac_pos) .+ 1e-8
+    center = SA[0.5, 0.5, 0.5]
+    wrapped_atom_frac_pos .+= center
+
+    # Remap fractional coordinates to atom positions
+    wrapped_atom_pos = get_real_lattice(atoms) * wrapped_atom_frac_pos
+
+    # Construct new atoms
+    recentered_atoms = periodic_system(
+        [Atom(atom, position = wrapped_atom_pos[:, iat]) for (iat, atom) in enumerate(atoms)],
+        cell_vectors(atoms),
+    )
+    return recentered_atoms
+end
+
+function get_real_lattice(atoms::AbstractSystem)
+    return stack(cell_vectors(atoms))
+end
+
+function get_recip_lattice(atoms::AbstractSystem)
+    real_lattice = get_real_lattice(atoms)
+    return get_recip_lattice(real_lattice)
+end
+
+# A^T B = 2πI
+# where both matrices contain basis vectors as columns
+function get_recip_lattice(real_lattice)
+    return transpose(ustrip.(real_lattice)) \ (2π * I(3))
+end
+
+# Strip units for matrix division because it doesn't always work with units.
+# Resulting fractional coordinates do not have units which is correct
+function get_frac_positions(atoms::AbstractSystem)
+    real_lattice = get_real_lattice(atoms)
+    real_positions = ustrip.(stack(position(atoms, :)))
+    return ustrip.(real_lattice) \ ustrip.(real_positions)
+end