Add symmetry

src/ExactDiagonalization.jl

@@ -7,6 +7,9 @@ using SparseArrays
 include("util.jl")
 include("HilbertSpace.jl")
 include("Operator.jl")
+
+include("Symmetry.jl")
+
 include("prettyprintln.jl")
 
 end

src/Symmetry.jl

@@ -1,9 +1,6 @@
-module Symmetry
 
 include("Symmetry/permutation.jl")
 include("Symmetry/group.jl")
 include("Symmetry/translation.jl")
 include("Symmetry/application.jl")
-
-
-end
+include("Symmetry/reduction.jl")

src/Symmetry/application.jl

@@ -1,3 +1,5 @@
+export apply_symmetry
+export is_invariant
 
 function apply_symmetry(hs::HilbertSpace, permutation ::Permutation, bitrep ::BR) where {BR}
   out = zero(BR)
@@ -7,3 +9,34 @@ function apply_symmetry(hs::HilbertSpace, permutation ::Permutation, bitrep ::BR
   return out
 end
 
+
+function apply_symmetry(permutation ::Permutation, op::NullOperator)
+  return op
+end
+
+
+function apply_symmetry(permutation ::Permutation, op::PureOperator{S, BR}) where {S,BR}
+  hs = op.hilbert_space
+  bm = apply_symmetry(op.hilbert_space, permutation, op.bitmask)
+  br = apply_symmetry(op.hilbert_space, permutation, op.bitrow)
+  bc = apply_symmetry(op.hilbert_space, permutation, op.bitcol)
+  am = op.amplitude
+  return PureOperator{S, BR}(hs, bm, br, bc, am)
+end
+
+
+function apply_symmetry(permutation ::Permutation, op::SumOperator{S, BR}) where {S, BR}
+  terms = collect(apply_symmetry(permutation, t) for t in op.terms)
+  return SumOperator{S, BR}(op.hilbert_space, terms)
+end
+
+
+function is_invariant(permutation ::Permutation, op::AbstractOperator)
+  return simplify(op - apply_symmetry(permutation, op)) == NullOperator()
+end
+
+function is_invariant(trans_group ::TranslationGroup, op::AbstractOperator)
+  return all(
+    is_invariant(g, op) for g in trans_group.generators
+  )
+end

src/Symmetry/perm.jl → src/Symmetry/permutation.jl

@@ -1,5 +1,10 @@
+export AbstractSymmetryOperation
+export Permutation
 
-struct Permutation
+
+abstract type AbstractSymmetryOperation end
+
+struct Permutation <: AbstractSymmetryOperation
   map ::Vector{Int}
   cycle_length ::Int
   function Permutation(perms ::AbstractVector{Int}; max_cycle=2048)

src/Symmetry/reduction.jl

@@ -0,0 +1,108 @@
+export ReducedHilbertSpaceRealization
+export symmetry_reduce
+export materialize
+
+struct ReducedHilbertSpaceRealization
+  parent_hilbert_space_realization ::HilbertSpaceRealization
+  translation_group ::TranslationGroup
+  basis_list ::Vector{UInt}
+  basis_lookup ::Dict{UInt, NamedTuple{(:index, :amplitude), Tuple{Int, ComplexF64}}}
+end
+
+function symmetry_reduce(hsr ::HilbertSpaceRealization,
+                trans_group ::TranslationGroup,
+                fractional_momentum ::AbstractVector{Rational})
+  ik = findfirst(collect(
+    trans_group.fractional_momenta[ik] == fractional_momentum
+    for ik in 1:length(trans_group.fractional_momenta) ))
+
+  isnothing(ik) && throw(ArgumentError("fractional momentum $(fractional_momentum) not an irrep of the translation group"))
+
+  # check if fractional momentum is compatible with translation group ?
+  #k = float.(fractional_momentum) .* 2π
+  phases = trans_group.character_table[ik, :]
+  #[ cis(dot(k, t)) for t in trans_group.translations]
+  reduced_basis_list = Set{UInt}()  
+  parent_amplitude = Dict()
+
+  for bvec in hsr.basis_list
+    if haskey(parent_amplitude, bvec)
+      continue
+    end
+
+    ψ = SparseState{ComplexF64, UInt}(hsr.hilbert_space)
+    identity_translations = Vector{Int}[]
+    for i in 1:length(trans_group.elements)
+      t = trans_group.translations[i]
+      g = trans_group.elements[i]
+      p = phases[i]
+
+      bvec_prime = apply_symmetry(hsr.hilbert_space, g, bvec)
+      ψ[bvec_prime] += p
+      if bvec_prime == bvec
+        push!(identity_translations, t)
+      end
+    end
+    if !is_compatible(fractional_momentum, identity_translations)
+      continue
+    end
+    clean!(ψ)
+    @assert !isempty(ψ)
+
+    normalize!(ψ)
+    push!(reduced_basis_list, bvec)
+
+    for (bvec_prime, amplitude) in ψ.components
+      parent_amplitude[bvec_prime] = (parent=bvec, amplitude=amplitude)
+    end
+  end
+  reduced_basis_list = sort(collect(reduced_basis_list))
+  reduced_basis_lookup = Dict(bvec => (index=ivec, amplitude=parent_amplitude[bvec].amplitude)
+                              for (ivec, bvec) in enumerate(reduced_basis_list))
+
+  for (bvec_prime, (bvec, amplitude)) in parent_amplitude
+    bvec_prime == bvec && continue
+    reduced_basis_lookup[bvec_prime] = (index=reduced_basis_lookup[bvec].index, amplitude=amplitude)
+  end
+  return ReducedHilbertSpaceRealization(hsr, trans_group, reduced_basis_list, reduced_basis_lookup)
+end
+
+
+function materialize(rhsr :: ReducedHilbertSpaceRealization,
+                     operator ::AbstractOperator;
+                     tol::Real=sqrt(eps(Float64)))
+  # TODO CHECK IF THe OPERATOR HAS TRANSLATION SYMMETRY
+  rows = Int[]
+  cols = Int[]
+  vals = ComplexF64[]
+
+  hs = rhsr.parent_hilbert_space_realization.hilbert_space
+  err = 0.0
+
+  for (irow, brow) in enumerate(rhsr.basis_list)
+    ampl_row = rhsr.basis_lookup[brow].amplitude
+    ψrow = SparseState{ComplexF64, UInt}(hs, brow=>1/ampl_row)
+    ψcol = SparseState{ComplexF64, UInt}(hs)
+    apply!(ψcol, ψrow, operator)
+    clean!(ψcol)
+
+    for (bcol, ampl) in ψcol.components
+      if ! haskey(rhsr.basis_lookup, bcol)
+        err += abs(ampl.^2)
+        continue
+      end
+
+      (icol, ampl_col) = rhsr.basis_lookup[bcol]
+      push!(rows, irow)
+      push!(cols, icol)
+      push!(vals, ampl * ampl_col)
+    end
+  end
+
+  if maximum(imag.(vals)) < tol
+    vals = real.(vals)
+  end
+
+  n = length(rhsr.basis_list)
+  return (sparse(rows, cols, vals, n, n), err)
+end

src/Symmetry/translation.jl

@@ -1,29 +1,46 @@
+export TranslationGroup
+export is_compatible
 
 struct TranslationGroup <: AbstractSymmetryGroup
   generators ::Vector{Permutation}
 
   translations ::Vector{Vector{Int}}
   elements ::Vector{Permutation}
   fractional_momenta ::Vector{Vector{Rational}}
-  #representations ::Vector{ Vector{ComplexF64} }
 
+  conjugacy_classes ::Vector{Int}  # conjugacy class of element
+  character_table ::Array{ComplexF64, 2}
+
   function TranslationGroup(generators::AbstractArray{Permutation})
-    for g1 in generators, g2 in generators
-      @assert g1 * g2 == g2 * g1
-    end
+    @assert all(g1 * g2 == g2 * g1 for g1 in generators, g2 in generators)
+
     shape = [g.cycle_length for g in generators]
     translations = vcat( collect( Iterators.product([0:g.cycle_length-1 for g in generators]...) )...)
     translations = [ [x...] for x in translations]
     elements = [prod(gen^d for (gen, d) in zip(generators, dist)) for (ig, dist) in enumerate(translations)]
+
     momentum(sub) = [x//d for (x, d) in zip(sub, shape)]
     fractional_momenta = [momentum(sub) for sub in translations]
-    #momentum(sub) = [x//d for (x, d) in zip(sub, shape)]
-    #momenta = [momentum(sub) for sub in translations]
-    #representations = collect(Iterators.product([0:g.cycle_length-1 for g in generators]...))
-    return new(generators, translations, elements, fractional_momenta)
+
+    conjugacy_classes = collect(1:length(elements))
+    character_table = ComplexF64[cis(dot(float.(kf) .* 2π, t))
+                                 for kf in fractional_momenta, t in translations]
+
+    return new(generators, translations, elements, fractional_momenta, conjugacy_classes, character_table)
   end
 end
 
+
+"""
+    is_compatible
+
+Check whether the fractional momentum ([0, 1)ᴺ) compatible with the identity translation.
+i.e. k¹ R¹ + k² R² + ... + kᴺ Rᴺ = 0 (mod 1)
+
+# Arguments
+- `fractional_momentum ::AbstractVector{Rational}` : k
+- `identity_translation ::AbstractVector{<:Integer}` : R
+"""
 function is_compatible(
     fractional_momentum ::AbstractVector{Rational},
     identity_translation ::AbstractVector{<:Integer}