From 43b5225fee28e189233bb608bd5d6f83b964b617 Mon Sep 17 00:00:00 2001 From: victor Date: Wed, 29 Jan 2025 10:54:05 +0100 Subject: [PATCH 1/3] add `ClusterExpansion` time evolution MPOs --- src/MPSKit.jl | 3 +- src/algorithms/timestep/clusterexpansion.jl | 148 ++++++++++++++++++++ 2 files changed, 150 insertions(+), 1 deletion(-) create mode 100644 src/algorithms/timestep/clusterexpansion.jl diff --git a/src/MPSKit.jl b/src/MPSKit.jl index 1fe569db2..1c11394bf 100644 --- a/src/MPSKit.jl +++ b/src/MPSKit.jl @@ -34,7 +34,7 @@ export VUMPS, VOMPS, DMRG, DMRG2, IDMRG, IDMRG2, GradientGrassmann export excitations export FiniteExcited, QuasiparticleAnsatz, ChepigaAnsatz, ChepigaAnsatz2 export time_evolve, timestep, timestep!, make_time_mpo -export TDVP, TDVP2, WI, WII, TaylorCluster +export TDVP, TDVP2, WI, WII, TaylorCluster, ClusterExpansion export changebonds, changebonds! export VUMPSSvdCut, OptimalExpand, SvdCut, RandExpand export propagator @@ -144,6 +144,7 @@ include("algorithms/changebonds/randexpand.jl") include("algorithms/timestep/tdvp.jl") include("algorithms/timestep/timeevmpo.jl") +include("algorithms/timestep/clusterexpansion.jl") include("algorithms/timestep/integrators.jl") include("algorithms/timestep/time_evolve.jl") diff --git a/src/algorithms/timestep/clusterexpansion.jl b/src/algorithms/timestep/clusterexpansion.jl new file mode 100644 index 000000000..2bd64cd59 --- /dev/null +++ b/src/algorithms/timestep/clusterexpansion.jl @@ -0,0 +1,148 @@ +struct ClusterExpansion + N::Int +end + +function make_time_mpo(H::MPOHamiltonian, dt::Number, alg::ClusterExpansion) + N = alg.N + lmax = N ÷ 2 # largest level + τ = -im * dt + # spaces + P = physicalspace(H)[1] + D = dim(physicalspace(H)[1]) # physical dimension + V = BlockTensorKit.oplus([ℂ^(D^2l) for l in 0:lmax]...) + + TT = tensormaptype(ComplexSpace, 2, 2, ComplexF64) + O = SparseBlockTensorMap{TT}(undef, V ⊗ P ← P ⊗ V) + + for n in 1:N + if n == 1 + O[1, 1, 1, 1] = _make_Hamterms(H, 1, τ) + elseif n == 2 + U, S, V = tsvd(_make_b(H, O, 2, τ), ((1, 2, 4), (3, 5, 6))) + O[1, 1, 1, 2] = permute(U * sqrt(S), ((1, 2), (3, 4))) + O[2, 1, 1, 1] = permute(sqrt(S) * V, ((1, 2), (3, 4))) + else + nl = n ÷ 2 # new level + jnl = nl + 1 # Julia indexing + apply_function = _make_apply_functions(O, n) + b = _make_b(H, O, n, τ) + Onew, info = lssolve(apply_function, b; verbosity=5) + + (info.converged) != 1 && @warn "Did not converge when constucting the $n cluster" + + if isodd(n) # linear problem + # assign the new level to O + O[jnl, 1, 1, jnl] = Onew + elseif iseven(n) # linear problem + svd + U, S, V = tsvd(Onew, (1, 2, 4), (3, 5, 6)) + # assign the new levels to O + O[jnl-1, 1, 1, jnl] = permute(U * sqrt(S), ((1, 2), (3, 4))) + O[jnl, 1, 1, jnl-1] = permute(sqrt(S) * V, ((1, 2), (3, 4))) + end + end + end + + # return O + return MPO(PeriodicVector([O])) +end + + +function _make_Hamterms(H::MPOHamiltonian, n::Int, τ::Number) + return add_util_leg(convert(TensorMap, scale(DenseMPO(open_boundary_conditions(H, n)), τ))) +end + +function _make_b(H::MPOHamiltonian, O::SparseBlockTensorMap, n::Int, τ::Number) + # be - bo + be = permute(exp(permute(_make_Hamterms(H, n, τ), Tuple(1:n+1), Tuple(vcat([2n + 2,], collect(n+2:2n+1))))), Tuple(1:n+1), Tuple(vcat(collect(n+3:2n+2), [n + 2,]))) + bo = _fully_contract_O(O, n) + return be - bo +end + +function _fully_contract_O(O::SparseBlockTensorMap, n::Int) + # Make projector onto the 0 subspace of the virtual level + Pl = zeros(ComplexF64, ℂ^1 ← space(O, 1)) + Pl[1] = 1 + + Pr = zeros(ComplexF64, space(O, 4)' ← ℂ^1) + Pr[1] = 1 + + O_inds = [[i, -(i + 1), -(i + 1 + n), i + 1] for i in 1:n] + + # Contract and permute to right structure + return permute(ncon([Pl, fill(O, n)..., Pr], [[-1, 1], O_inds..., [n + 1, -(2n + 2)]]), Tuple(1:n+1), Tuple(n+2:2n+2)) +end + +function make_envs(O::SparseBlockTensorMap, n::Int) + nt = (n - 1) ÷ 2 # amount of tensors in the environment + return _make_envs(O, nt) +end + +function _make_envs(O::SparseBlockTensorMap, nt::Int) + if nt == 1 + return O[1, 1, 1, 2], O[2, 1, 1, 1] + else + Tl, Tr = _make_envs(O, nt - 1) + Ol = O[nt, 1, 1, nt+1] + Or = O[nt+1, 1, 1, nt] + return _make_left_env(Tl, Ol), _make_right_env(Or, Tr) + end +end + +@generated function _make_left_env(Tl::AbstractTensorMap, Ol::AbstractTensorMap) + N₁ = numin(Tl) + N = numind(Tl) + t_out = tensorexpr(:new_e_l, -(1:N₁+1), -(N₁+2:N+2)) + t_left = tensorexpr(:Tl, -(1:N₁), ((-(N₁+2:N)...), 1)) + t_right = tensorexpr(:Ol, (1, -(N₁ + 1)), (-(N + 1), -(N + 2))) + return macroexpand(@__MODULE__, + :(return @tensor $t_out := $t_left * $t_right)) +end + +@generated function _make_right_env(Or::AbstractTensorMap, Tr::AbstractTensorMap) + N₁ = numin(Tr) + N = numind(Tr) + t_out = tensorexpr(:new_e_r, -(1:N₁+1), -(N₁+2:N+2)) + t_left = tensorexpr(:Or, -(1:2), (-(N₁ + 2), 1)) + t_right = tensorexpr(:Tr, (1, (-(3:N₁+1)...)), -(N₁+3:N+2)) + return macroexpand(@__MODULE__, + :(return @tensor $t_out := $t_left * $t_right)) +end + +function _make_apply_functions(O::SparseBlockTensorMap, n::Int) + nT = (n - 1) ÷ 2 + Al, Ar = make_envs(O, n) + + fun = if iseven(n) + let + function A(x::TensorMap, ::Val{false}) + Al_inds = vcat(-collect(1:nT+1), -collect(2nT+4:3nT+3), [1,]) + x_inds = [1, -(nT + 2), -(nT + 3), -(3nT + 4), -(3nT + 5), 2] + Ar_inds = vcat([2,], -collect(nT+4:2nT+3), -collect(3nT+6:4nT+6)) + return permute(ncon([Al, x, Ar], [Al_inds, x_inds, Ar_inds]), Tuple(1:2nT+3), Tuple((2nT+4):(4nT+6))) + end + + function A(b::TensorMap, ::Val{true}) + Al_inds = vcat(collect(3:nT+2), [-1, 1], collect(nT+3:2nT+2)) + b_inds = vcat([1,], collect(nT+3:2nT+2), [-2, -3], collect(3nT+3:4nT+2), collect(3:nT+2), [-4, -5], collect(2nT+3:3nT+2), [2,]) + Ar_inds = vcat(collect(2nT+3:3nT+2), [2, -6], collect(3nT+3:4nT+2)) + return permute(ncon([adjoint(Al), b, adjoint(Ar)], [Al_inds, b_inds, Ar_inds]), (1, 2, 3), (4, 5, 6)) + end + end + elseif isodd(n) + let + function A(x::TensorMap, ::Val{false}) + Al_inds = vcat(-collect(1:nT+1), -collect(2nT+3:3nT+2), [1,]) + x_inds = [1, -(nT + 2), -(3nT + 3), 2] + Ar_inds = vcat([2,], -collect(nT+3:2nT+2), -collect(3nT+4:4nT+4)) + return permute(ncon([Al, x, Ar], [Al_inds, x_inds, Ar_inds]), Tuple(1:2nT+2), Tuple((2nT+3):(4nT+4))) + end + function A(b::TensorMap, ::Val{true}) + Al_inds = vcat(collect(3:nT+2), [-1, 1], collect(nT+3:2nT+2)) + b_inds = vcat([1,], collect(nT+3:2nT+2), [-2,], collect(3nT+3:4nT+2), collect(3:nT+2), [-3,], collect(2nT+3:3nT+2), [2,]) + Ar_inds = vcat(collect(2nT+3:3nT+2), [2, -4], collect(3nT+3:4nT+2)) + return permute(ncon([adjoint(Al), b, adjoint(Ar)], [Al_inds, b_inds, Ar_inds]), (1, 2), (3, 4)) + end + end + end + return fun +end \ No newline at end of file From 4b8af2c03f6fb52607e7359d2b9f2ae90086ae7b Mon Sep 17 00:00:00 2001 From: victor Date: Tue, 18 Feb 2025 14:17:53 +0100 Subject: [PATCH 2/3] format --- src/algorithms/timestep/clusterexpansion.jl | 76 +++++++++++++-------- 1 file changed, 46 insertions(+), 30 deletions(-) diff --git a/src/algorithms/timestep/clusterexpansion.jl b/src/algorithms/timestep/clusterexpansion.jl index 2bd64cd59..2260e20c8 100644 --- a/src/algorithms/timestep/clusterexpansion.jl +++ b/src/algorithms/timestep/clusterexpansion.jl @@ -28,7 +28,8 @@ function make_time_mpo(H::MPOHamiltonian, dt::Number, alg::ClusterExpansion) b = _make_b(H, O, n, τ) Onew, info = lssolve(apply_function, b; verbosity=5) - (info.converged) != 1 && @warn "Did not converge when constucting the $n cluster" + (info.converged) != 1 && + @warn "Did not converge when constucting the $n cluster" if isodd(n) # linear problem # assign the new level to O @@ -36,8 +37,8 @@ function make_time_mpo(H::MPOHamiltonian, dt::Number, alg::ClusterExpansion) elseif iseven(n) # linear problem + svd U, S, V = tsvd(Onew, (1, 2, 4), (3, 5, 6)) # assign the new levels to O - O[jnl-1, 1, 1, jnl] = permute(U * sqrt(S), ((1, 2), (3, 4))) - O[jnl, 1, 1, jnl-1] = permute(sqrt(S) * V, ((1, 2), (3, 4))) + O[jnl - 1, 1, 1, jnl] = permute(U * sqrt(S), ((1, 2), (3, 4))) + O[jnl, 1, 1, jnl - 1] = permute(sqrt(S) * V, ((1, 2), (3, 4))) end end end @@ -46,14 +47,16 @@ function make_time_mpo(H::MPOHamiltonian, dt::Number, alg::ClusterExpansion) return MPO(PeriodicVector([O])) end - function _make_Hamterms(H::MPOHamiltonian, n::Int, τ::Number) - return add_util_leg(convert(TensorMap, scale(DenseMPO(open_boundary_conditions(H, n)), τ))) + return add_util_leg(convert(TensorMap, + scale(DenseMPO(open_boundary_conditions(H, n)), τ))) end function _make_b(H::MPOHamiltonian, O::SparseBlockTensorMap, n::Int, τ::Number) # be - bo - be = permute(exp(permute(_make_Hamterms(H, n, τ), Tuple(1:n+1), Tuple(vcat([2n + 2,], collect(n+2:2n+1))))), Tuple(1:n+1), Tuple(vcat(collect(n+3:2n+2), [n + 2,]))) + be = permute(exp(permute(_make_Hamterms(H, n, τ), Tuple(1:(n + 1)), + Tuple(vcat([2n + 2], collect((n + 2):(2n + 1)))))), + Tuple(1:(n + 1)), Tuple(vcat(collect((n + 3):(2n + 2)), [n + 2]))) bo = _fully_contract_O(O, n) return be - bo end @@ -69,7 +72,8 @@ function _fully_contract_O(O::SparseBlockTensorMap, n::Int) O_inds = [[i, -(i + 1), -(i + 1 + n), i + 1] for i in 1:n] # Contract and permute to right structure - return permute(ncon([Pl, fill(O, n)..., Pr], [[-1, 1], O_inds..., [n + 1, -(2n + 2)]]), Tuple(1:n+1), Tuple(n+2:2n+2)) + return permute(ncon([Pl, fill(O, n)..., Pr], [[-1, 1], O_inds..., [n + 1, -(2n + 2)]]), + Tuple(1:(n + 1)), Tuple((n + 2):(2n + 2))) end function make_envs(O::SparseBlockTensorMap, n::Int) @@ -82,8 +86,8 @@ function _make_envs(O::SparseBlockTensorMap, nt::Int) return O[1, 1, 1, 2], O[2, 1, 1, 1] else Tl, Tr = _make_envs(O, nt - 1) - Ol = O[nt, 1, 1, nt+1] - Or = O[nt+1, 1, 1, nt] + Ol = O[nt, 1, 1, nt + 1] + Or = O[nt + 1, 1, 1, nt] return _make_left_env(Tl, Ol), _make_right_env(Or, Tr) end end @@ -91,21 +95,21 @@ end @generated function _make_left_env(Tl::AbstractTensorMap, Ol::AbstractTensorMap) N₁ = numin(Tl) N = numind(Tl) - t_out = tensorexpr(:new_e_l, -(1:N₁+1), -(N₁+2:N+2)) - t_left = tensorexpr(:Tl, -(1:N₁), ((-(N₁+2:N)...), 1)) + t_out = tensorexpr(:new_e_l, -(1:(N₁ + 1)), -((N₁ + 2):(N + 2))) + t_left = tensorexpr(:Tl, -(1:N₁), ((-((N₁ + 2):N)...), 1)) t_right = tensorexpr(:Ol, (1, -(N₁ + 1)), (-(N + 1), -(N + 2))) return macroexpand(@__MODULE__, - :(return @tensor $t_out := $t_left * $t_right)) + :(return @tensor $t_out := $t_left * $t_right)) end @generated function _make_right_env(Or::AbstractTensorMap, Tr::AbstractTensorMap) N₁ = numin(Tr) N = numind(Tr) - t_out = tensorexpr(:new_e_r, -(1:N₁+1), -(N₁+2:N+2)) + t_out = tensorexpr(:new_e_r, -(1:(N₁ + 1)), -((N₁ + 2):(N + 2))) t_left = tensorexpr(:Or, -(1:2), (-(N₁ + 2), 1)) - t_right = tensorexpr(:Tr, (1, (-(3:N₁+1)...)), -(N₁+3:N+2)) + t_right = tensorexpr(:Tr, (1, (-(3:(N₁ + 1))...)), -((N₁ + 3):(N + 2))) return macroexpand(@__MODULE__, - :(return @tensor $t_out := $t_left * $t_right)) + :(return @tensor $t_out := $t_left * $t_right)) end function _make_apply_functions(O::SparseBlockTensorMap, n::Int) @@ -115,34 +119,46 @@ function _make_apply_functions(O::SparseBlockTensorMap, n::Int) fun = if iseven(n) let function A(x::TensorMap, ::Val{false}) - Al_inds = vcat(-collect(1:nT+1), -collect(2nT+4:3nT+3), [1,]) + Al_inds = vcat(-collect(1:(nT + 1)), -collect((2nT + 4):(3nT + 3)), [1]) x_inds = [1, -(nT + 2), -(nT + 3), -(3nT + 4), -(3nT + 5), 2] - Ar_inds = vcat([2,], -collect(nT+4:2nT+3), -collect(3nT+6:4nT+6)) - return permute(ncon([Al, x, Ar], [Al_inds, x_inds, Ar_inds]), Tuple(1:2nT+3), Tuple((2nT+4):(4nT+6))) + Ar_inds = vcat([2], -collect((nT + 4):(2nT + 3)), + -collect((3nT + 6):(4nT + 6))) + return permute(ncon([Al, x, Ar], [Al_inds, x_inds, Ar_inds]), + Tuple(1:(2nT + 3)), Tuple((2nT + 4):(4nT + 6))) end function A(b::TensorMap, ::Val{true}) - Al_inds = vcat(collect(3:nT+2), [-1, 1], collect(nT+3:2nT+2)) - b_inds = vcat([1,], collect(nT+3:2nT+2), [-2, -3], collect(3nT+3:4nT+2), collect(3:nT+2), [-4, -5], collect(2nT+3:3nT+2), [2,]) - Ar_inds = vcat(collect(2nT+3:3nT+2), [2, -6], collect(3nT+3:4nT+2)) - return permute(ncon([adjoint(Al), b, adjoint(Ar)], [Al_inds, b_inds, Ar_inds]), (1, 2, 3), (4, 5, 6)) + Al_inds = vcat(collect(3:(nT + 2)), [-1, 1], collect((nT + 3):(2nT + 2))) + b_inds = vcat([1], collect((nT + 3):(2nT + 2)), [-2, -3], + collect((3nT + 3):(4nT + 2)), collect(3:(nT + 2)), [-4, -5], + collect((2nT + 3):(3nT + 2)), [2]) + Ar_inds = vcat(collect((2nT + 3):(3nT + 2)), [2, -6], + collect((3nT + 3):(4nT + 2))) + return permute(ncon([adjoint(Al), b, adjoint(Ar)], + [Al_inds, b_inds, Ar_inds]), (1, 2, 3), (4, 5, 6)) end end elseif isodd(n) let function A(x::TensorMap, ::Val{false}) - Al_inds = vcat(-collect(1:nT+1), -collect(2nT+3:3nT+2), [1,]) + Al_inds = vcat(-collect(1:(nT + 1)), -collect((2nT + 3):(3nT + 2)), [1]) x_inds = [1, -(nT + 2), -(3nT + 3), 2] - Ar_inds = vcat([2,], -collect(nT+3:2nT+2), -collect(3nT+4:4nT+4)) - return permute(ncon([Al, x, Ar], [Al_inds, x_inds, Ar_inds]), Tuple(1:2nT+2), Tuple((2nT+3):(4nT+4))) + Ar_inds = vcat([2], -collect((nT + 3):(2nT + 2)), + -collect((3nT + 4):(4nT + 4))) + return permute(ncon([Al, x, Ar], [Al_inds, x_inds, Ar_inds]), + Tuple(1:(2nT + 2)), Tuple((2nT + 3):(4nT + 4))) end function A(b::TensorMap, ::Val{true}) - Al_inds = vcat(collect(3:nT+2), [-1, 1], collect(nT+3:2nT+2)) - b_inds = vcat([1,], collect(nT+3:2nT+2), [-2,], collect(3nT+3:4nT+2), collect(3:nT+2), [-3,], collect(2nT+3:3nT+2), [2,]) - Ar_inds = vcat(collect(2nT+3:3nT+2), [2, -4], collect(3nT+3:4nT+2)) - return permute(ncon([adjoint(Al), b, adjoint(Ar)], [Al_inds, b_inds, Ar_inds]), (1, 2), (3, 4)) + Al_inds = vcat(collect(3:(nT + 2)), [-1, 1], collect((nT + 3):(2nT + 2))) + b_inds = vcat([1], collect((nT + 3):(2nT + 2)), [-2], + collect((3nT + 3):(4nT + 2)), collect(3:(nT + 2)), [-3], + collect((2nT + 3):(3nT + 2)), [2]) + Ar_inds = vcat(collect((2nT + 3):(3nT + 2)), [2, -4], + collect((3nT + 3):(4nT + 2))) + return permute(ncon([adjoint(Al), b, adjoint(Ar)], + [Al_inds, b_inds, Ar_inds]), (1, 2), (3, 4)) end end end return fun -end \ No newline at end of file +end From 179687725dd77d3d6869c023978d4ea41103a7bd Mon Sep 17 00:00:00 2001 From: victor Date: Tue, 24 Jun 2025 15:52:02 +0200 Subject: [PATCH 3/3] apply suggestions --- src/algorithms/timestep/clusterexpansion.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/algorithms/timestep/clusterexpansion.jl b/src/algorithms/timestep/clusterexpansion.jl index 2260e20c8..23181ef2b 100644 --- a/src/algorithms/timestep/clusterexpansion.jl +++ b/src/algorithms/timestep/clusterexpansion.jl @@ -7,8 +7,8 @@ function make_time_mpo(H::MPOHamiltonian, dt::Number, alg::ClusterExpansion) lmax = N ÷ 2 # largest level τ = -im * dt # spaces - P = physicalspace(H)[1] - D = dim(physicalspace(H)[1]) # physical dimension + P = physicalspace(H, 1) + D = dim(P) # physical dimension V = BlockTensorKit.oplus([ℂ^(D^2l) for l in 0:lmax]...) TT = tensormaptype(ComplexSpace, 2, 2, ComplexF64)