From 9392a9f540e3cf6b9b8204a4776412119f3b33ec Mon Sep 17 00:00:00 2001 From: Kyungmin Lee Date: Tue, 1 Oct 2019 11:01:38 +0900 Subject: [PATCH] Add symmetry --- src/ExactDiagonalization.jl | 3 + src/Symmetry.jl | 5 +- src/Symmetry/application.jl | 33 +++++++ src/Symmetry/{perm.jl => permutation.jl} | 7 +- src/Symmetry/reduction.jl | 108 +++++++++++++++++++++++ src/Symmetry/translation.jl | 33 +++++-- 6 files changed, 176 insertions(+), 13 deletions(-) rename src/Symmetry/{perm.jl => permutation.jl} (90%) create mode 100644 src/Symmetry/reduction.jl diff --git a/src/ExactDiagonalization.jl b/src/ExactDiagonalization.jl index f14f07d..be5a948 100644 --- a/src/ExactDiagonalization.jl +++ b/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 diff --git a/src/Symmetry.jl b/src/Symmetry.jl index b3752bf..8e18c87 100644 --- a/src/Symmetry.jl +++ b/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 \ No newline at end of file +include("Symmetry/reduction.jl") \ No newline at end of file diff --git a/src/Symmetry/application.jl b/src/Symmetry/application.jl index 682960f..b410db8 100644 --- a/src/Symmetry/application.jl +++ b/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 \ No newline at end of file diff --git a/src/Symmetry/perm.jl b/src/Symmetry/permutation.jl similarity index 90% rename from src/Symmetry/perm.jl rename to src/Symmetry/permutation.jl index 4f9155f..50f8491 100644 --- a/src/Symmetry/perm.jl +++ b/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) diff --git a/src/Symmetry/reduction.jl b/src/Symmetry/reduction.jl new file mode 100644 index 0000000..1f8b8fa --- /dev/null +++ b/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 diff --git a/src/Symmetry/translation.jl b/src/Symmetry/translation.jl index 693b158..68aa235 100644 --- a/src/Symmetry/translation.jl +++ b/src/Symmetry/translation.jl @@ -1,3 +1,5 @@ +export TranslationGroup +export is_compatible struct TranslationGroup <: AbstractSymmetryGroup generators ::Vector{Permutation} @@ -5,25 +7,40 @@ struct TranslationGroup <: AbstractSymmetryGroup 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}