Generalize infsiteinds to fractional QN densities

src/infinitecanonicalmps.jl

@@ -12,11 +12,23 @@ function Base.:/(qnval::ITensors.QNVal, n::Int)
   return setval(qnval, Int(div_val))
 end
 
+function Base.:*(qnval::ITensors.QNVal, n::Int)
+  div_val = ITensors.val(qnval) * n
+  if !isinteger(div_val)
+    error("Multiplying $qnval by $n, the resulting QN value is not an integer")
+  end
+  return setval(qnval, Int(div_val))
+end
+
 # TODO: Move to ITensors.jl
 function Base.:/(qn::QN, n::Int)
   return QN(map(qnval -> qnval / n, qn.data))
 end
 
+function Base.:*(qn::QN, n::Int)
+  return QN(map(qnval -> qnval * n, qn.data))
+end
+
 # of Index (Tuple, Vector, ITensor, etc.)
 indtype(i::Index) = typeof(i)
 indtype(T::Type{<:Index}) = T
@@ -48,11 +60,82 @@ function shift_flux(i::Index, flux_density::QN)
   return ITensors.setspace(i, shift_flux(space(i), flux_density))
 end
 
+function multiply_flux(qn::ITensors.QNVal, flux_factor::Int64)
+  if qn.modulus == 0 #is undefined, so remains undefined
+    return (qn.name, qn.val, qn.modulus) #qn
+  elseif qn.modulus == 1 || qn.modulus == -1 #U(1) symmetry, should remain a U(1) symmetry
+    return (qn.name, qn.val * flux_factor, qn.modulus) #qn * flux_factor
+  else #To get the same algebra, multiplying the modulus by flux_factor works
+    return (qn.name, qn.val * flux_factor, qn.modulus * flux_factor)
+  end
+end
+
+function multiply_flux(qn::QN, flux_factor::Int64)
+  flux_factor == 1 && return qn
+  return QN(multiply_flux.(filter(x -> x.modulus != 0, qn.data), flux_factor)...)
+end
+
+function multiply_flux(qn::QN, flux_factor::Dict{ITensors.SmallString,Int64})
+  return QN(
+    [multiply_flux(q, flux_factor[q.name]) for q in filter(x -> x.modulus != 0, qn.data)]...
+  )
+end
+
+function multiply_flux(
+  qnblock::Pair{QN,Int}, flux_factor::Union{Int64,Dict{ITensors.SmallString,Int64}}
+)
+  flux_factor == 1 && return qnblock
+  return (multiply_flux(ITensors.qn(qnblock), flux_factor) => ITensors.blockdim(qnblock))
+end
+
+function multiply_flux(
+  space::Vector{Pair{QN,Int}}, flux_factor::Union{Int64,Dict{ITensors.SmallString,Int64}}
+)
+  flux_factor == 1 && return space
+  return map(qnblock -> multiply_flux(qnblock, flux_factor), space)
+end
+
+function multiply_flux(i::Index, flux_factor::Union{Int64,Dict{ITensors.SmallString,Int64}})
+  flux_factor == 1 && return i
+  return ITensors.setspace(i, multiply_flux(space(i), flux_factor))
+end
+
+#= Alternative version where we do not try to optimize the multiplier for each QN
 function shift_flux_to_zero(s::Vector{<:Index}, flux::QN)
   if iszero(flux)
     return s
   end
+  #For simplicity, we are not introducing a factor per QN
   n = length(s)
+  multiplier = 1
+  for qn in flux
+    if qn.val == 0
+      continue
+    end
+    multiplier = lcm(lcm(qn.val, n)÷qn.val, multiplier)
+  end
+  s = map(sₙ -> multiply_flux(sₙ, multiplier), s)
+  flux = multiply_flux(flux, multiplier)
+  flux_density = flux / n
+  return map(sₙ -> shift_flux(sₙ, flux_density), s)
+end
+=#
+
+function shift_flux_to_zero(s::Vector{<:Index}, flux::QN)
+  if iszero(flux)
+    return s
+  end
+  #For simplicity, we are not introducing a factor per QN
+  n = length(s)
+  multipliers = Dict{ITensors.SmallString,Int64}() #multipliers is assigned using the names
+  for qn in flux.data
+    if qn.modulus == 0
+      continue
+    end
+    multipliers[qn.name] = qn.val != 0 ? lcm(qn.val, n) ÷ qn.val : 1
+  end
+  s = map(sₙ -> multiply_flux(sₙ, multipliers), s)
+  flux = multiply_flux(flux, multipliers)
   flux_density = flux / n
   return map(sₙ -> shift_flux(sₙ, flux_density), s)
 end

src/subspace_expansion.jl

@@ -18,7 +18,7 @@ function generate_twobody_nullspace(
 
   range_H = nrange(ψ, H[1])
 
-  @assert range_H > 1 "Not defined for purely local Hamiltonians"
+  @assert range_H > 1 "Subspace expansion for InfiniteSum{MPO} is not defined for purely local Hamiltonians"
 
   if range_H > 2
     ψᴴ = dag(ψ)
@@ -181,7 +181,7 @@ function subspace_expansion(
   @assert dˡ == dʳ
   if dˡ ≥ maxdim
     println(
-      "Current bond dimension at bond $b is $dˡ while desired maximum dimension is $maxdim, skipping bond dimension increase",
+      "Current bond dimension at bond $b is $dˡ while desired maximum dimension is $maxdim, skipping bond dimension increase at $b",
     )
     flush(stdout)
     flush(stderr)
@@ -199,7 +199,7 @@ function subspace_expansion(
   #Added due to crash during testing
   if norm(ψHN2.tensor) < 1e-12
     println(
-      "Impossible to do a subspace expansion, probably due to conservation constraints"
+      "The two-site subspace expansion produced a zero-norm expansion at $b. This is likely due to the long-range nature of the QN conserving Hamiltonian.",
     )
     flush(stdout)
     flush(stderr)

test/test_vumpsmpo.jl

@@ -9,7 +9,7 @@ function expect_three_site(ψ::MPS, h::ITensor, n::Int)
 end
 
 #Ref time is 21.6s with negligible compilation time
-@time @testset "vumpsmpo_ising" begin
+@testset "vumpsmpo_ising" begin
   Random.seed!(1234)
 
   model = Model("ising")

test/test_vumpsmpo_fqhe.jl

@@ -21,7 +21,7 @@ end
 
 #Currently, VUMPS cannot give the right result as the subspace expansion is too small
 #This is meant to test the generalized translation rules
-@time @testset "vumpsmpo_fqhe" begin
+@testset "vumpsmpo_fqhe" begin
   Random.seed!(1234)
   function initstate(n)
     if mod(n, 3) == 1