From 4ff1e6ddb31f24b0a3126e41e52905e59abe2ff2 Mon Sep 17 00:00:00 2001 From: CompatHelper Julia Date: Tue, 24 Jan 2023 01:14:12 +0000 Subject: [PATCH 01/17] CompatHelper: bump compat for KrylovKit to 0.6, (keep existing compat) --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 7d2104d..c77ea83 100644 --- a/Project.toml +++ b/Project.toml @@ -22,7 +22,7 @@ HDF5 = "0.15, 0.16" ITensors = "0.3.18" Infinities = "0.1" IterTools = "1" -KrylovKit = "0.5" +KrylovKit = "0.5, 0.6" OffsetArrays = "1" QuadGK = "2" Requires = "1" From 2df126c5bfb3f58f6673aa1026b6698ebe84d3d0 Mon Sep 17 00:00:00 2001 From: mtfishman Date: Mon, 20 Feb 2023 01:10:45 +0000 Subject: [PATCH 02/17] Format .jl files --- src/ITensorInfiniteMPS.jl | 3 ++- test/test_vumpsmpo.jl | 3 ++- test/test_vumpsmpo_fqhe.jl | 3 ++- 3 files changed, 6 insertions(+), 3 deletions(-) diff --git a/src/ITensorInfiniteMPS.jl b/src/ITensorInfiniteMPS.jl index 996cfbb..a5fabee 100644 --- a/src/ITensorInfiniteMPS.jl +++ b/src/ITensorInfiniteMPS.jl @@ -87,7 +87,8 @@ export Cell, function __init__() # This is used for debugging using visualizations - @require ITensorsVisualization = "f2aed53d-2f32-47c3-a7b9-1ee253853786" @eval using ITensorsVisualization + @require ITensorsVisualization = "f2aed53d-2f32-47c3-a7b9-1ee253853786" + @eval using ITensorsVisualization end end diff --git a/test/test_vumpsmpo.jl b/test/test_vumpsmpo.jl index 98b1185..97714a8 100644 --- a/test/test_vumpsmpo.jl +++ b/test/test_vumpsmpo.jl @@ -9,7 +9,8 @@ function expect_three_site(ψ::MPS, h::ITensor, n::Int) end #Ref time is 21.6s with negligible compilation time -@time @testset "vumpsmpo_ising" begin +@time +@testset "vumpsmpo_ising" begin Random.seed!(1234) model = Model("ising") diff --git a/test/test_vumpsmpo_fqhe.jl b/test/test_vumpsmpo_fqhe.jl index 76a9d22..c6a2c92 100644 --- a/test/test_vumpsmpo_fqhe.jl +++ b/test/test_vumpsmpo_fqhe.jl @@ -21,7 +21,8 @@ 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 +@time +@testset "vumpsmpo_fqhe" begin Random.seed!(1234) function initstate(n) if mod(n, 3) == 1 From 92d48c78d1730fb74524e54df6d58dfd543f1028 Mon Sep 17 00:00:00 2001 From: LHerviou <98113223+LHerviou@users.noreply.github.com> Date: Mon, 20 Feb 2023 17:16:15 +0100 Subject: [PATCH 03/17] Fixing infsiteinds to accept "fractional" fluxes --- src/infinitecanonicalmps.jl | 19 +++++++++++++++++-- 1 file changed, 17 insertions(+), 2 deletions(-) diff --git a/src/infinitecanonicalmps.jl b/src/infinitecanonicalmps.jl index 2831ba0..ec582de 100644 --- a/src/infinitecanonicalmps.jl +++ b/src/infinitecanonicalmps.jl @@ -48,13 +48,28 @@ function shift_flux(i::Index, flux_density::QN) return ITensors.setspace(i, shift_flux(space(i), flux_density)) end +function multiply_flux(qnblock::Pair{QN,Int}, flux_factor::Int64) + return ((ITensors.qn(qnblock)) => ITensors.blockdim(qnblock)) +end +function multiply_flux(space::Vector{Pair{QN,Int}}, flux_factor::Int64) + return map(qnblock -> multiply_flux(qnblock, flux_factor), space) +end +function multiply_flux(i::Index, flux_factor::Int64) + return ITensors.setspace(i, multiply_flux(space(i), flux_factor)) +end + function shift_flux_to_zero(s::Vector{<:Index}, flux::QN) if iszero(flux) return s end n = length(s) - flux_density = flux / n - return map(sₙ -> shift_flux(sₙ, flux_density), s) + try + flux_density = flux / n + return map(sₙ -> shift_flux(sₙ, flux_density), s) + catch e + s = map(sₙ -> multiply_flux(sₙ, n), s) + return map(sₙ -> shift_flux(sₙ, flux), s) + end end function infsiteinds( From 1ef65892165804475625ad066c8073c174da3757 Mon Sep 17 00:00:00 2001 From: LHerviou <98113223+LHerviou@users.noreply.github.com> Date: Mon, 20 Feb 2023 17:39:10 +0100 Subject: [PATCH 04/17] Update ITensorInfiniteMPS.jl --- src/ITensorInfiniteMPS.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/ITensorInfiniteMPS.jl b/src/ITensorInfiniteMPS.jl index a5fabee..3afa5e9 100644 --- a/src/ITensorInfiniteMPS.jl +++ b/src/ITensorInfiniteMPS.jl @@ -87,8 +87,7 @@ export Cell, function __init__() # This is used for debugging using visualizations - @require ITensorsVisualization = "f2aed53d-2f32-47c3-a7b9-1ee253853786" - @eval using ITensorsVisualization + @require ITensorsVisualization="f2aed53d-2f32-47c3-a7b9-1ee253853786" @eval using ITensorsVisualization end end From c99718d098559c1eccf3003bfe3696d60b9b75f6 Mon Sep 17 00:00:00 2001 From: LHerviou Date: Mon, 20 Feb 2023 17:44:39 +0100 Subject: [PATCH 05/17] Fixing formatting --- src/ITensorInfiniteMPS.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ITensorInfiniteMPS.jl b/src/ITensorInfiniteMPS.jl index 3afa5e9..996cfbb 100644 --- a/src/ITensorInfiniteMPS.jl +++ b/src/ITensorInfiniteMPS.jl @@ -87,7 +87,7 @@ export Cell, function __init__() # This is used for debugging using visualizations - @require ITensorsVisualization="f2aed53d-2f32-47c3-a7b9-1ee253853786" @eval using ITensorsVisualization + @require ITensorsVisualization = "f2aed53d-2f32-47c3-a7b9-1ee253853786" @eval using ITensorsVisualization end end From 5a53cbdb072d8b2ba75633787563987994dcfa89 Mon Sep 17 00:00:00 2001 From: LHerviou Date: Mon, 20 Feb 2023 18:19:57 +0100 Subject: [PATCH 06/17] Forgot to define multiply --- src/infinitecanonicalmps.jl | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/src/infinitecanonicalmps.jl b/src/infinitecanonicalmps.jl index ec582de..3e63f4f 100644 --- a/src/infinitecanonicalmps.jl +++ b/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 @@ -49,7 +61,7 @@ function shift_flux(i::Index, flux_density::QN) end function multiply_flux(qnblock::Pair{QN,Int}, flux_factor::Int64) - return ((ITensors.qn(qnblock)) => ITensors.blockdim(qnblock)) + return ((ITensors.qn(qnblock)*flux_factor) => ITensors.blockdim(qnblock)) end function multiply_flux(space::Vector{Pair{QN,Int}}, flux_factor::Int64) return map(qnblock -> multiply_flux(qnblock, flux_factor), space) From cf767cfa1c4cbf9f64c999b7046b326ae5c2fe94 Mon Sep 17 00:00:00 2001 From: mtfishman Date: Tue, 21 Feb 2023 01:11:42 +0000 Subject: [PATCH 07/17] Format .jl files --- src/infinitecanonicalmps.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/infinitecanonicalmps.jl b/src/infinitecanonicalmps.jl index 3e63f4f..70df054 100644 --- a/src/infinitecanonicalmps.jl +++ b/src/infinitecanonicalmps.jl @@ -61,7 +61,7 @@ function shift_flux(i::Index, flux_density::QN) end function multiply_flux(qnblock::Pair{QN,Int}, flux_factor::Int64) - return ((ITensors.qn(qnblock)*flux_factor) => ITensors.blockdim(qnblock)) + return ((ITensors.qn(qnblock) * flux_factor) => ITensors.blockdim(qnblock)) end function multiply_flux(space::Vector{Pair{QN,Int}}, flux_factor::Int64) return map(qnblock -> multiply_flux(qnblock, flux_factor), space) From 86e20645f2d17f6490a3fc3c663264caaa06ef7c Mon Sep 17 00:00:00 2001 From: LHerviou Date: Mon, 20 Mar 2023 15:43:18 +0100 Subject: [PATCH 08/17] Polished treatment of qns --- src/infinitecanonicalmps.jl | 76 ++++++++++++++++++++++++++++++++----- 1 file changed, 66 insertions(+), 10 deletions(-) diff --git a/src/infinitecanonicalmps.jl b/src/infinitecanonicalmps.jl index 01c25ba..b65d919 100644 --- a/src/infinitecanonicalmps.jl +++ b/src/infinitecanonicalmps.jl @@ -60,28 +60,84 @@ function shift_flux(i::Index, flux_density::QN) return ITensors.setspace(i, shift_flux(space(i), flux_density)) end -function multiply_flux(qnblock::Pair{QN,Int}, flux_factor::Int64) - return ((ITensors.qn(qnblock) * flux_factor) => ITensors.blockdim(qnblock)) +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(space::Vector{Pair{QN,Int}}, flux_factor::Int64) + +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::Int64) + +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) - try - flux_density = flux / n - return map(sₙ -> shift_flux(sₙ, flux_density), s) - catch e - s = map(sₙ -> multiply_flux(sₙ, n), s) - return map(sₙ -> shift_flux(sₙ, flux), 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 function infsiteinds( From aa290feeab433dac6ee5964b18125c204b959131 Mon Sep 17 00:00:00 2001 From: mtfishman Date: Tue, 21 Mar 2023 01:03:26 +0000 Subject: [PATCH 09/17] Format .jl files --- src/infinitecanonicalmps.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/infinitecanonicalmps.jl b/src/infinitecanonicalmps.jl index 7d1a8ec..01c25ba 100644 --- a/src/infinitecanonicalmps.jl +++ b/src/infinitecanonicalmps.jl @@ -61,7 +61,7 @@ function shift_flux(i::Index, flux_density::QN) end function multiply_flux(qnblock::Pair{QN,Int}, flux_factor::Int64) - return ((ITensors.qn(qnblock)*flux_factor) => ITensors.blockdim(qnblock)) + return ((ITensors.qn(qnblock) * flux_factor) => ITensors.blockdim(qnblock)) end function multiply_flux(space::Vector{Pair{QN,Int}}, flux_factor::Int64) return map(qnblock -> multiply_flux(qnblock, flux_factor), space) From 73f1f18aa286e2a831ab888219431c12f115b111 Mon Sep 17 00:00:00 2001 From: LHerviou Date: Tue, 21 Mar 2023 09:45:05 +0100 Subject: [PATCH 10/17] Fixing tests --- test/test_vumpsmpo.jl | 1 - test/test_vumpsmpo_fqhe.jl | 1 - 2 files changed, 2 deletions(-) diff --git a/test/test_vumpsmpo.jl b/test/test_vumpsmpo.jl index 97714a8..ea2d705 100644 --- a/test/test_vumpsmpo.jl +++ b/test/test_vumpsmpo.jl @@ -9,7 +9,6 @@ function expect_three_site(ψ::MPS, h::ITensor, n::Int) end #Ref time is 21.6s with negligible compilation time -@time @testset "vumpsmpo_ising" begin Random.seed!(1234) diff --git a/test/test_vumpsmpo_fqhe.jl b/test/test_vumpsmpo_fqhe.jl index c6a2c92..98a156f 100644 --- a/test/test_vumpsmpo_fqhe.jl +++ b/test/test_vumpsmpo_fqhe.jl @@ -21,7 +21,6 @@ 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 Random.seed!(1234) function initstate(n) From 79739c4a5c6f55989dcc3b78dfcc3eadf7104c72 Mon Sep 17 00:00:00 2001 From: LHerviou Date: Tue, 21 Mar 2023 12:38:16 +0100 Subject: [PATCH 11/17] Improving some error messages --- src/subspace_expansion.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/subspace_expansion.jl b/src/subspace_expansion.jl index c0bd39a..5d2eb20 100644 --- a/src/subspace_expansion.jl +++ b/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) From 695e0ced2988d63a374f1f858ab38883cb16c95d Mon Sep 17 00:00:00 2001 From: LHerviou Date: Wed, 22 Mar 2023 14:27:39 +0100 Subject: [PATCH 12/17] Taking into account Mat comment --- src/infinitecanonicalmps.jl | 76 ++++++++++++------------------------- 1 file changed, 24 insertions(+), 52 deletions(-) diff --git a/src/infinitecanonicalmps.jl b/src/infinitecanonicalmps.jl index b65d919..81739f5 100644 --- a/src/infinitecanonicalmps.jl +++ b/src/infinitecanonicalmps.jl @@ -13,11 +13,8 @@ function Base.:/(qnval::ITensors.QNVal, n::Int) 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)) + prod_val = ITensors.val(qnval) * n + return setval(qnval, Int(prod_val)) end # TODO: Move to ITensors.jl @@ -60,66 +57,41 @@ 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 +function scale_flux(qn::ITensors.QNVal, flux_factor::Int) + flux_factor == 1 && return qn + if qn.modulus == 0 #should not be called directly, but is a better safety + return qn #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 + return 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) + return ITensors.QNVal(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)...) +function scale_flux(qn::ITensors.QNVal, flux_factor::Dict{ITensors.SmallString,Int}) + qn.name == ITensors.SmallString() && return qn #Ensure that we do not touch empty QNs + return scale_flux(qn, flux_factor[qn.name]) 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)]... - ) +function scale_flux(qn::QN, flux_factor::Union{Int,Dict{ITensors.SmallString,Int}}) + return QN(map(qnval -> scale_flux(qnval, flux_factor), qn.data)) end -function multiply_flux( - qnblock::Pair{QN,Int}, flux_factor::Union{Int64,Dict{ITensors.SmallString,Int64}} +function scale_flux( + qnblock::Pair{QN,Int}, flux_factor::Union{Int,Dict{ITensors.SmallString,Int}} ) - flux_factor == 1 && return qnblock - return (multiply_flux(ITensors.qn(qnblock), flux_factor) => ITensors.blockdim(qnblock)) + return (scale_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}} +function scale_flux( + space::Vector{Pair{QN,Int}}, flux_factor::Union{Int,Dict{ITensors.SmallString,Int}} ) - flux_factor == 1 && return space - return map(qnblock -> multiply_flux(qnblock, flux_factor), space) + return map(qnblock -> scale_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) +function scale_flux(i::Index, flux_factor::Union{Int,Dict{ITensors.SmallString,Int}}) + return ITensors.setspace(i, scale_flux(space(i), flux_factor)) end -=# function shift_flux_to_zero(s::Vector{<:Index}, flux::QN) if iszero(flux) @@ -127,15 +99,15 @@ function shift_flux_to_zero(s::Vector{<:Index}, flux::QN) 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 + multipliers = Dict{ITensors.SmallString,Int}() #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) + s = map(sₙ -> scale_flux(sₙ, multipliers), s) + flux = scale_flux(flux, multipliers) flux_density = flux / n return map(sₙ -> shift_flux(sₙ, flux_density), s) end From 60f979977d73ad8e9ca49749a921189a82131bcd Mon Sep 17 00:00:00 2001 From: LHerviou <98113223+LHerviou@users.noreply.github.com> Date: Wed, 22 Mar 2023 15:28:23 +0100 Subject: [PATCH 13/17] Update infinitecanonicalmps.jl --- src/infinitecanonicalmps.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/infinitecanonicalmps.jl b/src/infinitecanonicalmps.jl index 81739f5..fd2ca54 100644 --- a/src/infinitecanonicalmps.jl +++ b/src/infinitecanonicalmps.jl @@ -104,7 +104,7 @@ function shift_flux_to_zero(s::Vector{<:Index}, flux::QN) if qn.modulus == 0 continue end - multipliers[qn.name] = qn.val != 0 ? lcm(qn.val, n) ÷ qn.val : 1 + multipliers[qn.name] = qn.val != 0 ? abs(lcm(qn.val, n) ÷ qn.val) : 1 #It is a bit more readable to keep multipliers positive. end s = map(sₙ -> scale_flux(sₙ, multipliers), s) flux = scale_flux(flux, multipliers) From eb1dacacbc075c0e0bcd80f8a4aa7ba39dbf8f9d Mon Sep 17 00:00:00 2001 From: LHerviou Date: Wed, 22 Mar 2023 15:54:12 +0100 Subject: [PATCH 14/17] Further simplification of the code --- src/infinitecanonicalmps.jl | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/src/infinitecanonicalmps.jl b/src/infinitecanonicalmps.jl index fd2ca54..49d01f1 100644 --- a/src/infinitecanonicalmps.jl +++ b/src/infinitecanonicalmps.jl @@ -58,10 +58,8 @@ function shift_flux(i::Index, flux_density::QN) end function scale_flux(qn::ITensors.QNVal, flux_factor::Int) - flux_factor == 1 && return qn - if qn.modulus == 0 #should not be called directly, but is a better safety - return qn #qn - elseif qn.modulus == 1 || qn.modulus == -1 #U(1) symmetry, should remain a U(1) symmetry + flux_factor == 1 && return qn #This test could occur earlier + if -1 <= qn.modulus <= 1 #U(1) symmetry, should remain a U(1) symmetry, if qn.modulus = 0, it is a dummy term and the code does nothing return qn * flux_factor else #To get the same algebra, multiplying the modulus by flux_factor works return ITensors.QNVal(qn.name, qn.val * flux_factor, qn.modulus * flux_factor) @@ -69,7 +67,7 @@ function scale_flux(qn::ITensors.QNVal, flux_factor::Int) end function scale_flux(qn::ITensors.QNVal, flux_factor::Dict{ITensors.SmallString,Int}) - qn.name == ITensors.SmallString() && return qn #Ensure that we do not touch empty QNs + !haskey(flux_factor, qn.name) && return qn #Ensure that we do not touch empty QNs or QNs that we do not need to rescale because the shift is 0 return scale_flux(qn, flux_factor[qn.name]) end @@ -97,14 +95,14 @@ 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 + #We are introducing a factor per QN n = length(s) multipliers = Dict{ITensors.SmallString,Int}() #multipliers is assigned using the names for qn in flux.data - if qn.modulus == 0 + if qn.val == 0 #We do not need to multiply in this case. It also takes into account the empty QNs which have val and modulus 0 continue end - multipliers[qn.name] = qn.val != 0 ? abs(lcm(qn.val, n) ÷ qn.val) : 1 #It is a bit more readable to keep multipliers positive. + multipliers[qn.name] = abs(lcm(qn.val, n) ÷ qn.val) #It is a bit more readable to keep multipliers positive. end s = map(sₙ -> scale_flux(sₙ, multipliers), s) flux = scale_flux(flux, multipliers) From 4801dbf9c6e08aa7b510cedd8cb9d12949cdd0e3 Mon Sep 17 00:00:00 2001 From: LHerviou Date: Mon, 27 Mar 2023 15:41:26 +0200 Subject: [PATCH 15/17] Better auto qns for Zn symmetries --- src/infinitecanonicalmps.jl | 44 ++++++++++++++++++++++++++++++---- test/test_unitcell_qns.jl | 47 +++++++++++++++++++++++++++++++++++++ 2 files changed, 86 insertions(+), 5 deletions(-) create mode 100644 test/test_unitcell_qns.jl diff --git a/src/infinitecanonicalmps.jl b/src/infinitecanonicalmps.jl index 49d01f1..4e43e41 100644 --- a/src/infinitecanonicalmps.jl +++ b/src/infinitecanonicalmps.jl @@ -5,11 +5,31 @@ end # TODO: Move to ITensors.jl function Base.:/(qnval::ITensors.QNVal, n::Int) - div_val = ITensors.val(qnval) / n - if !isinteger(div_val) - error("Dividing $qnval by $n, the resulting QN value is not an integer") + if abs(qnval.modulus) <= 1 + div_val = ITensors.val(qnval) / n + if !isinteger(div_val) + error("Dividing $qnval by $n, the resulting QN value is not an integer") + end + return setval(qnval, Int(div_val)) + else + if mod(qnval.val, gcd(qnval.modulus, n)) != 0 + error("Dividing $qnval by $n, no solution to the Chinese remainder theorem") + end + #We perform a brute force sieving solution -> should be ok for all reasonnable n + sol = qnval.val + for x in 1:n + if mod(sol, n) != 0 + sol += qnval.modulus + else + break + end + x == n && error("Failed to find solution") #Bezut identity would be a better solver + end + modulus = lcm(qnval.modulus, n) + sol = mod(sol, modulus) + sol = sol < abs(sol - modulus) ? sol : sol - modulus + return setval(qnval, Int(sol / n)) end - return setval(qnval, Int(div_val)) end function Base.:*(qnval::ITensors.QNVal, n::Int) @@ -102,7 +122,21 @@ function shift_flux_to_zero(s::Vector{<:Index}, flux::QN) if qn.val == 0 #We do not need to multiply in this case. It also takes into account the empty QNs which have val and modulus 0 continue end - multipliers[qn.name] = abs(lcm(qn.val, n) ÷ qn.val) #It is a bit more readable to keep multipliers positive. + multipliers[qn.name] = abs(lcm(qn.val, n) ÷ qn.val)#default solution with positive multiplier + if abs(qn.modulus) > 1 + #We find a multiplier which guarantees a solution of the Chinese theorem + #For now, brute force: + multipliers[qn.name] = abs(lcm(qn.val, n) ÷ qn.val) + for x in 1:(multipliers[qn.name] - 1) + if mod(n, x) != 0 + continue + end + if mod(qn.val * x, gcd(n, x * qn.modulus)) == 0 + multipliers[qn.name] = x + break + end + end + end end s = map(sₙ -> scale_flux(sₙ, multipliers), s) flux = scale_flux(flux, multipliers) diff --git a/test/test_unitcell_qns.jl b/test/test_unitcell_qns.jl new file mode 100644 index 0000000..da58084 --- /dev/null +++ b/test/test_unitcell_qns.jl @@ -0,0 +1,47 @@ +using ITensors +using ITensorInfiniteMPS +using Test +using Random + +@testset "unitcell_qns" begin + function initstate(n) + if n % 2 == 0 + if (n / 2) % 2 == 1 + return "Dn" + else + return "Up" + end + else + if (n + 1) / 2 <= 19 + if ((n + 1) / 2) % 2 == 1 + return "Up" + else + return "Dn" + end + else + return "Emp" + end + end + end + conserve_qns = true + N = 62 + s = infsiteinds(n -> isodd(n) ? "tJ" : "S=1/2", N; conserve_qns, initstate) + ψ = InfMPS(s, initstate) + @show N + @show s + + function initstate(n) + if mod(n, 4) == 1 + return 2 + else + return 1 + end + end + + for N in 1:6 + s = infsiteinds("S=1/2", N; conserve_szparity=true, initstate) + ψ = InfMPS(s, initstate) + @show N + @show s + end +end From 1a6038bf00bafce353e9b95ebe1eff2eef3f10c8 Mon Sep 17 00:00:00 2001 From: LHerviou Date: Tue, 28 Mar 2023 16:10:27 +0200 Subject: [PATCH 16/17] Comments and test --- src/infinitecanonicalmps.jl | 28 +++++++++++++++++++++------- test/test_unitcell_qns.jl | 4 +++- 2 files changed, 24 insertions(+), 8 deletions(-) diff --git a/src/infinitecanonicalmps.jl b/src/infinitecanonicalmps.jl index 4e43e41..a32da24 100644 --- a/src/infinitecanonicalmps.jl +++ b/src/infinitecanonicalmps.jl @@ -15,6 +15,8 @@ function Base.:/(qnval::ITensors.QNVal, n::Int) if mod(qnval.val, gcd(qnval.modulus, n)) != 0 error("Dividing $qnval by $n, no solution to the Chinese remainder theorem") end + #We look for the inverse of n in the equation n x = qnval.val mod qn.modulus. + #The Chinese remainder theorem guarantees it exists #We perform a brute force sieving solution -> should be ok for all reasonnable n sol = qnval.val for x in 1:n @@ -23,7 +25,7 @@ function Base.:/(qnval::ITensors.QNVal, n::Int) else break end - x == n && error("Failed to find solution") #Bezut identity would be a better solver + x == n && error("Failed to find solution") #Bezut identity would be a better solver here end modulus = lcm(qnval.modulus, n) sol = mod(sol, modulus) @@ -124,15 +126,27 @@ function shift_flux_to_zero(s::Vector{<:Index}, flux::QN) end multipliers[qn.name] = abs(lcm(qn.val, n) ÷ qn.val)#default solution with positive multiplier if abs(qn.modulus) > 1 - #We find a multiplier which guarantees a solution of the Chinese theorem - #For now, brute force: + #= + for Zp symmetries, we use the Chinese remainder theorem to find an optimal solution. + We try to find y such that y = 0 mod n and y = qn.val mod qn.modulus. Then we just need to find x such that y = n x. + + This system admits a solution iff qn.val = 0 mod( gcd(qn.modulus, n)), which is not necessarily true. + In order to avoid this problem, we introduce a global multiplier b: qn.val -> b qn.val, qn.modulus -> b qn.modulus. + If b = n, we always have n qn.val = 0 mod( gcd(n qn.modulus, n) = n ) + We can do better: I try to find the smallest b such that: + b qn.val = 0 mod gcd(n, qn.modulus b). + + Case when things changed: n = 3, initstate = +--, flux = 1 mod 2. + Previous version: 1 is not divisible by 3, so we mutiply everything by 3. Gives QN ( 4 mod 6, 1 mod 6) + Current version: gcd(2, 3) = 1, so a solution actually exists: 3 x 1 = 1 mod 2 => work with QN (0 mod 2, 1 mod 2) + =# multipliers[qn.name] = abs(lcm(qn.val, n) ÷ qn.val) - for x in 1:(multipliers[qn.name] - 1) - if mod(n, x) != 0 + for b in 1:(multipliers[qn.name] - 1) + if mod(n, b) != 0 #It is easy to show that the optimal b divides n continue end - if mod(qn.val * x, gcd(n, x * qn.modulus)) == 0 - multipliers[qn.name] = x + if mod(qn.val * b, gcd(n, b * qn.modulus)) == 0 + multipliers[qn.name] = b break end end diff --git a/test/test_unitcell_qns.jl b/test/test_unitcell_qns.jl index da58084..a2caa26 100644 --- a/test/test_unitcell_qns.jl +++ b/test/test_unitcell_qns.jl @@ -40,8 +40,10 @@ using Random for N in 1:6 s = infsiteinds("S=1/2", N; conserve_szparity=true, initstate) - ψ = InfMPS(s, initstate) @show N @show s + @test iszero(flux(s)) + ψ = InfMPS(s, initstate) + @test iszero(flux(ψ)) end end From 7234efa319089d93af44eb713bd5a551df9e8b9d Mon Sep 17 00:00:00 2001 From: LHerviou Date: Tue, 28 Mar 2023 16:14:05 +0200 Subject: [PATCH 17/17] Debug --- test/test_unitcell_qns.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/test/test_unitcell_qns.jl b/test/test_unitcell_qns.jl index a2caa26..3102fa0 100644 --- a/test/test_unitcell_qns.jl +++ b/test/test_unitcell_qns.jl @@ -42,8 +42,7 @@ using Random s = infsiteinds("S=1/2", N; conserve_szparity=true, initstate) @show N @show s - @test iszero(flux(s)) ψ = InfMPS(s, initstate) - @test iszero(flux(ψ)) + @test iszero(flux(ψ.AL)) end end