diff --git a/src/algorithms/toolbox.jl b/src/algorithms/toolbox.jl index 95582b66f..466338541 100644 --- a/src/algorithms/toolbox.jl +++ b/src/algorithms/toolbox.jl @@ -36,8 +36,8 @@ function infinite_temperature_density_matrix(H::MPOHamiltonian) end """ - calc_galerkin(above, operator, below, envs) - calc_galerkin(pos, above, operator, below, envs) + calc_galerkin(below, operator, above, envs) + calc_galerkin(pos, below, operator, above, envs) Calculate the Galerkin error, which is the error between the solution of the original problem, and the solution of the problem projected on the tangent space. Concretely, this is the overlap of the current state with the single-site derivative, projected onto the nullspace of the current state: @@ -46,25 +46,25 @@ Concretely, this is the overlap of the current state with the single-site deriva \\epsilon = |VL * (VL' * \\frac{above}{\\partial AC_{pos}})| ``` """ -function calc_galerkin(pos::Int, above::Union{InfiniteMPS,FiniteMPS,WindowMPS}, operator, - below, envs) - AC´ = ∂∂AC(pos, above, operator, envs) * above.AC[pos] +function calc_galerkin(pos::Int, below::Union{InfiniteMPS,FiniteMPS,WindowMPS}, operator, + above, envs) + AC´ = ∂∂AC(pos, below, operator, envs) * above.AC[pos] normalize!(AC´) out = add!(AC´, below.AL[pos] * below.AL[pos]' * AC´, -1) return norm(out) end -function calc_galerkin(pos::CartesianIndex{2}, above::MultilineMPS, operator::MultilineMPO, - below::MultilineMPS, envs::MultilineEnvironments) +function calc_galerkin(pos::CartesianIndex{2}, below::MultilineMPS, operator::MultilineMPO, + above::MultilineMPS, envs::MultilineEnvironments) row, col = pos.I - return calc_galerkin(col, above[row], operator[row], below[row + 1], envs[row]) + return calc_galerkin(col, below[row + 1], operator[row], above[row], envs[row]) end -function calc_galerkin(above::Union{InfiniteMPS,FiniteMPS,WindowMPS}, operator, - below, envs) - return maximum(pos -> calc_galerkin(pos, above, operator, below, envs), 1:length(above)) +function calc_galerkin(below::Union{InfiniteMPS,FiniteMPS,WindowMPS}, operator, + above, envs) + return maximum(pos -> calc_galerkin(pos, below, operator, above, envs), 1:length(above)) end -function calc_galerkin(above::MultilineMPS, operator::MultilineMPO, below::MultilineMPS, +function calc_galerkin(below::MultilineMPS, operator::MultilineMPO, above::MultilineMPS, envs::MultilineEnvironments) - return maximum(pos -> calc_galerkin(pos, above, operator, below, envs), + return maximum(pos -> calc_galerkin(pos, below, operator, above, envs), CartesianIndices(size(above))) end diff --git a/src/environments/abstract_envs.jl b/src/environments/abstract_envs.jl index 01daac8c9..0efc0e1ce 100644 --- a/src/environments/abstract_envs.jl +++ b/src/environments/abstract_envs.jl @@ -64,7 +64,7 @@ end # Environment algorithms # ---------------------- """ - environment_alg(above, operator, below; kwargs...) + environment_alg(below, operator, above; kwargs...) Determine an appropriate algorithm for computing the environments, based on the given `kwargs...`. """ @@ -76,7 +76,7 @@ function environment_alg(::Union{InfiniteMPS,MultilineMPS}, eager=true) return Arnoldi(; tol, maxiter, krylovdim, verbosity, eager) end -function environment_alg(above, ::InfiniteMPOHamiltonian, below; +function environment_alg(below, ::InfiniteMPOHamiltonian, above; tol=Defaults.tol, maxiter=Defaults.maxiter, krylovdim=Defaults.krylovdim, verbosity=Defaults.VERBOSE_NONE) return GMRES(; tol, maxiter, krylovdim, verbosity) diff --git a/src/environments/infinite_envs.jl b/src/environments/infinite_envs.jl index d6ab383e4..62d74a1b2 100644 --- a/src/environments/infinite_envs.jl +++ b/src/environments/infinite_envs.jl @@ -18,19 +18,19 @@ Base.length(envs::InfiniteEnvironments) = length(envs.GLs) leftenv(envs::InfiniteEnvironments, site::Int, state) = envs.GLs[site] rightenv(envs::InfiniteEnvironments, site::Int, state) = envs.GRs[site] -function environments(above::InfiniteMPS, +function environments(below::InfiniteMPS, operator::Union{InfiniteMPO,InfiniteMPOHamiltonian}, - below::InfiniteMPS=above; + above::InfiniteMPS=below; kwargs...) - GLs, GRs = initialize_environments(above, operator, below) + GLs, GRs = initialize_environments(below, operator, above) envs = InfiniteEnvironments(GLs, GRs) - return recalculate!(envs, above, operator, below; kwargs...) + return recalculate!(envs, below, operator, above; kwargs...) end -function issamespace(envs::InfiniteEnvironments, above::InfiniteMPS, +function issamespace(envs::InfiniteEnvironments, below::InfiniteMPS, operator::Union{InfiniteMPO,InfiniteMPOHamiltonian}, - below::InfiniteMPS) - L = check_length(above, operator, below) + above::InfiniteMPS) + L = check_length(below, operator, above) for i in 1:L space(envs.GLs[i]) == (left_virtualspace(below, i) ⊗ left_virtualspace(operator, i)' ← @@ -42,39 +42,39 @@ function issamespace(envs::InfiniteEnvironments, above::InfiniteMPS, return true end -function recalculate!(envs::InfiniteEnvironments, above::InfiniteMPS, +function recalculate!(envs::InfiniteEnvironments, below::InfiniteMPS, operator::Union{InfiniteMPO,InfiniteMPOHamiltonian}, - below::InfiniteMPS=above; kwargs...) - if !issamespace(envs, above, operator, below) + above::InfiniteMPS=below; kwargs...) + if !issamespace(envs, below, operator, above) # TODO: in-place initialization? - GLs, GRs = initialize_environments(above, operator, below) + GLs, GRs = initialize_environments(below, operator, above) copy!(envs.GLs, GLs) copy!(envs.GRs, GRs) end - alg = environment_alg(above, operator, below; kwargs...) + alg = environment_alg(below, operator, above; kwargs...) @sync begin - @spawn compute_leftenvs!(envs, above, operator, below, alg) - @spawn compute_rightenvs!(envs, above, operator, below, alg) + @spawn compute_leftenvs!(envs, below, operator, above, alg) + @spawn compute_rightenvs!(envs, below, operator, above, alg) end - normalize!(envs, above, operator, below) + normalize!(envs, below, operator, above) return envs end # InfiniteMPO environments # ------------------------ -function initialize_environments(above::InfiniteMPS, operator::InfiniteMPO, - below::InfiniteMPS=above) - L = check_length(above, operator, below) +function initialize_environments(below::InfiniteMPS, operator::InfiniteMPO, + above::InfiniteMPS=below) + L = check_length(below, operator, above) GLs = PeriodicVector([randomize!(allocate_GL(below, operator, above, i)) for i in 1:L]) GRs = PeriodicVector([randomize!(allocate_GR(below, operator, above, i)) for i in 1:L]) return GLs, GRs end -function compute_leftenvs!(envs::InfiniteEnvironments, above::InfiniteMPS, - operator::InfiniteMPO, below::InfiniteMPS, alg) +function compute_leftenvs!(envs::InfiniteEnvironments, below::InfiniteMPS, + operator::InfiniteMPO, above::InfiniteMPS, alg) # compute eigenvector T = TransferMatrix(above.AL, operator, below.AL) λ, envs.GLs[1] = fixedpoint(flip(T), envs.GLs[1], :LM, alg) @@ -86,8 +86,8 @@ function compute_leftenvs!(envs::InfiniteEnvironments, above::InfiniteMPS, return λ, envs end -function compute_rightenvs!(envs::InfiniteEnvironments, above::InfiniteMPS, - operator::InfiniteMPO, below::InfiniteMPS, alg) +function compute_rightenvs!(envs::InfiniteEnvironments, below::InfiniteMPS, + operator::InfiniteMPO, above::InfiniteMPS, alg) # compute eigenvector T = TransferMatrix(above.AR, operator, below.AR) λ, envs.GRs[end] = fixedpoint(T, envs.GRs[end], :LM, alg) @@ -99,8 +99,8 @@ function compute_rightenvs!(envs::InfiniteEnvironments, above::InfiniteMPS, return λ, envs end -function TensorKit.normalize!(envs::InfiniteEnvironments, above::InfiniteMPS, - operator::InfiniteMPO, below::InfiniteMPS) +function TensorKit.normalize!(envs::InfiniteEnvironments, below::InfiniteMPS, + operator::InfiniteMPO, above::InfiniteMPS) for i in 1:length(operator) λ = dot(below.C[i], MPO_∂∂C(envs.GLs[i + 1], envs.GRs[i]) * above.C[i]) scale!(envs.GLs[i + 1], inv(λ)) @@ -110,8 +110,8 @@ end # InfiniteMPOHamiltonian environments # ----------------------------------- -function initialize_environments(above::InfiniteMPS, operator::InfiniteMPOHamiltonian, - below::InfiniteMPS=above) +function initialize_environments(below::InfiniteMPS, operator::InfiniteMPOHamiltonian, + above::InfiniteMPS=below) L = check_length(above, operator, below) GLs = PeriodicVector([allocate_GL(below, operator, above, i) for i in 1:L]) GRs = PeriodicVector([allocate_GR(below, operator, above, i) for i in 1:L]) @@ -139,9 +139,9 @@ function initialize_environments(above::InfiniteMPS, operator::InfiniteMPOHamilt return GLs, GRs end -function compute_leftenvs!(envs::InfiniteEnvironments, above::InfiniteMPS, - operator::InfiniteMPOHamiltonian, below::InfiniteMPS, alg) - L = check_length(above, operator, below) +function compute_leftenvs!(envs::InfiniteEnvironments, below::InfiniteMPS, + operator::InfiniteMPOHamiltonian, above::InfiniteMPS, alg) + L = check_length(below, above, operator) GLs = envs.GLs vsize = length(first(GLs)) @@ -156,12 +156,12 @@ function compute_leftenvs!(envs::InfiniteEnvironments, above::InfiniteMPS, # fill_data!(leftutil, one) # @plansor GL[1][1][-1 -2; -3] = ρ_left[-1; -3] * leftutil[-2] - (L > 1) && left_cyclethrough!(1, GLs, above, operator, below) + (L > 1) && left_cyclethrough!(1, GLs, below, operator, above) for i in 2:vsize prev = copy(GLs[1][i]) zerovector!(GLs[1][i]) - left_cyclethrough!(i, GLs, above, operator, below) + left_cyclethrough!(i, GLs, below, operator, above) if isidentitylevel(operator, i) # identity matrices; do the hacky renormalization T = regularize(TransferMatrix(above.AL, below.AL), ρ_left, ρ_right) @@ -169,7 +169,7 @@ function compute_leftenvs!(envs::InfiniteEnvironments, above::InfiniteMPS, convhist.converged == 0 && @warn "GL$i failed to converge: normres = $(convhist.normres)" - (L > 1) && left_cyclethrough!(i, GLs, above, operator, below) + (L > 1) && left_cyclethrough!(i, GLs, below, operator, above) # go through the unitcell, again subtracting fixpoints for site in 1:L @@ -186,15 +186,15 @@ function compute_leftenvs!(envs::InfiniteEnvironments, above::InfiniteMPS, convhist.converged == 0 && @warn "GL$i failed to converge: normres = $(convhist.normres)" end - (L > 1) && left_cyclethrough!(i, GLs, above, operator, below) + (L > 1) && left_cyclethrough!(i, GLs, below, operator, above) end end return GLs end -function left_cyclethrough!(index::Int, GL, above::InfiniteMPS, H::InfiniteMPOHamiltonian, - below::InfiniteMPS=above) +function left_cyclethrough!(index::Int, GL, below::InfiniteMPS, H::InfiniteMPOHamiltonian, + above::InfiniteMPS=below) # TODO: efficient transfer matrix slicing for large unitcells leftinds = 1:index for site in eachindex(GL) @@ -205,8 +205,8 @@ function left_cyclethrough!(index::Int, GL, above::InfiniteMPS, H::InfiniteMPOHa return GL end -function compute_rightenvs!(envs::InfiniteEnvironments, above::InfiniteMPS, - operator::InfiniteMPOHamiltonian, below::InfiniteMPS, alg) +function compute_rightenvs!(envs::InfiniteEnvironments, below::InfiniteMPS, + operator::InfiniteMPOHamiltonian, above::InfiniteMPS, alg) L = check_length(above, operator, below) GRs = envs.GRs vsize = length(last(GRs)) @@ -222,12 +222,12 @@ function compute_rightenvs!(envs::InfiniteEnvironments, above::InfiniteMPS, # fill_data!(rightutil, one) # @plansor GR[end][end][-1 -2; -3] = r_RR(state)[-1; -3] * rightutil[-2] - (L > 1) && right_cyclethrough!(vsize, GRs, above, operator, below) # populate other sites + (L > 1) && right_cyclethrough!(vsize, GRs, below, operator, above) # populate other sites for i in (vsize - 1):-1:1 prev = copy(GRs[end][i]) zerovector!(GRs[end][i]) - right_cyclethrough!(i, GRs, above, operator, below) + right_cyclethrough!(i, GRs, below, operator, above) if isidentitylevel(operator, i) # identity matrices; do the hacky renormalization # subtract fixpoints @@ -236,7 +236,7 @@ function compute_rightenvs!(envs::InfiniteEnvironments, above::InfiniteMPS, convhist.converged == 0 && @warn "GR$i failed to converge: normres = $(convhist.normres)" - L > 1 && right_cyclethrough!(i, GRs, above, operator, below) + L > 1 && right_cyclethrough!(i, GRs, below, operator, above) # go through the unitcell, again subtracting fixpoints for site in 1:L @@ -253,16 +253,16 @@ function compute_rightenvs!(envs::InfiniteEnvironments, above::InfiniteMPS, @warn "GR$i failed to converge: normres = $(convhist.normres)" end - (L > 1) && right_cyclethrough!(i, GRs, above, operator, below) + (L > 1) && right_cyclethrough!(i, GRs, below, operator, above) end end return GRs end -function right_cyclethrough!(index::Int, GR, above::InfiniteMPS, +function right_cyclethrough!(index::Int, GR, below::InfiniteMPS, operator::InfiniteMPOHamiltonian, - below::InfiniteMPS) + above::InfiniteMPS=below) # TODO: efficient transfer matrix slicing for large unitcells for site in reverse(eachindex(GR)) rightinds = index:length(GR[site]) @@ -274,21 +274,21 @@ function right_cyclethrough!(index::Int, GR, above::InfiniteMPS, end # no normalization necessary -- for consistant interface -function TensorKit.normalize!(envs::InfiniteEnvironments, above::InfiniteMPS, - operator::InfiniteMPOHamiltonian, below::InfiniteMPS) +function TensorKit.normalize!(envs::InfiniteEnvironments, below::InfiniteMPS, + operator::InfiniteMPOHamiltonian, above::InfiniteMPS) return envs end # Transfer operations # ------------------- -function transfer_leftenv!(envs::InfiniteEnvironments, above, operator, below, site::Int) +function transfer_leftenv!(envs::InfiniteEnvironments, below, operator, above, site::Int) T = TransferMatrix(above.AL[site - 1], operator[site - 1], below.AL[site - 1]) envs.GLs[site] = envs.GLs[site - 1] * T return envs end -function transfer_rightenv!(envs::InfiniteEnvironments, above, operator, below, site::Int) +function transfer_rightenv!(envs::InfiniteEnvironments, below, operator, above, site::Int) T = TransferMatrix(above.AR[site + 1], operator[site + 1], below.AR[site + 1]) envs.GRs[site] = T * envs.GRs[site + 1] return envs diff --git a/src/environments/multiline_envs.jl b/src/environments/multiline_envs.jl index de84fef0e..5d19666ce 100644 --- a/src/environments/multiline_envs.jl +++ b/src/environments/multiline_envs.jl @@ -1,38 +1,38 @@ const MultilineEnvironments{E<:AbstractMPSEnvironments} = Multiline{E} -function environments(above::MultilineMPS, operator::MultilineMPO, - below::MultilineMPS=above; kwargs...) +function environments(below::MultilineMPS, operator::MultilineMPO, + above::MultilineMPS=below; kwargs...) (rows = size(above, 1)) == size(operator, 1) == size(below, 1) || throw(ArgumentError("Incompatible sizes")) envs = map(1:rows) do row - return environments(above[row], operator[row], below[row + 1]; kwargs...) + return environments(below[row + 1], operator[row], above[row]; kwargs...) end return Multiline(PeriodicVector(envs)) end -function recalculate!(envs::MultilineEnvironments, above::MultilineMPS, - operator::MultilineMPO, below::MultilineMPS=above; kwargs...) +function recalculate!(envs::MultilineEnvironments, below::MultilineMPS, + operator::MultilineMPO, above::MultilineMPS=below; kwargs...) (rows = size(above, 1)) == size(operator, 1) == size(below, 1) || throw(ArgumentError("Incompatible sizes")) @threads for row in 1:rows - recalculate!(envs[row], above[row], operator[row], below[row + 1]; kwargs...) + recalculate!(envs[row], below[row + 1], operator[row], above[row]; kwargs...) end return envs end function recalculate!(envs::MultilineEnvironments, below, (operator, above)::Tuple; kwargs...) - return recalculate!(envs, above, operator, below; kwargs...) + return recalculate!(envs, below, operator, above; kwargs...) end -function TensorKit.normalize!(envs::MultilineEnvironments, above, operator, below) - for row in 1:size(above, 1) - normalize!(envs[row], above[row], operator[row], below[row + 1]) +function TensorKit.normalize!(envs::MultilineEnvironments, below, operator, above) + for row in 1:size(below, 1) + normalize!(envs[row], below[row + 1], operator[row], above[row]) end return envs end function TensorKit.normalize!(envs::MultilineEnvironments, below, (operator, above)) for row in 1:size(above, 1) - normalize!(envs[row], above[row], operator[row], below[row + 1]) + normalize!(envs[row], below[row + 1], operator[row], above[row]) end return envs end @@ -44,22 +44,22 @@ function rightenv(envs::MultilineEnvironments, col::Int, state) return rightenv.(parent(envs), col, parent(state)) end -function transfer_leftenv!(envs::MultilineEnvironments, above, operator, below, site::Int) +function transfer_leftenv!(envs::MultilineEnvironments, below, operator, above, site::Int) for row in 1:size(above, 1) - transfer_leftenv!(envs[row], above[row], operator[row], below[row + 1], site) + transfer_leftenv!(envs[row], below[row + 1], operator[row], above[row], site) end return envs end function transfer_leftenv!(envs::MultilineEnvironments, below, (O, above)::Tuple, site::Int) - return transfer_leftenv!(envs, above, O, below, site) + return transfer_leftenv!(envs, below, O, above, site) end -function transfer_rightenv!(envs::MultilineEnvironments, above, operator, below, site::Int) +function transfer_rightenv!(envs::MultilineEnvironments, below, operator, above, site::Int) for row in 1:size(above, 1) - transfer_rightenv!(envs[row], above[row], operator[row], below[row + 1], site) + transfer_rightenv!(envs[row], below[row + 1], operator[row], above[row], site) end return envs end function transfer_rightenv!(envs::MultilineEnvironments, below, (O, above)::Tuple, site::Int) - return transfer_rightenv!(envs, above, O, below, site) + return transfer_rightenv!(envs, below, O, above, site) end diff --git a/src/environments/multiple_envs.jl b/src/environments/multiple_envs.jl index 97a6f70f0..c1d806f5b 100644 --- a/src/environments/multiple_envs.jl +++ b/src/environments/multiple_envs.jl @@ -28,10 +28,10 @@ end """ Recalculate in-place each sub-env in MultipleEnvironments """ -function recalculate!(envs::MultipleEnvironments, above, operator::LazySum, below=above; +function recalculate!(envs::MultipleEnvironments, below, operator::LazySum, above=below; kwargs...) for (subenvs, subO) in zip(envs.envs, operator) - recalculate!(subenvs, above, subO, below; kwargs...) + recalculate!(subenvs, below, subO, above; kwargs...) end return envs end @@ -46,17 +46,17 @@ function Base.getproperty(envs::MultipleEnvironments, prop::Symbol) end function transfer_rightenv!(envs::MultipleEnvironments{<:InfiniteEnvironments}, - above, operator, below, pos::Int) + below, operator, above, pos::Int) for (subH, subenv) in zip(operator, envs.envs) - transfer_rightenv!(subenv, above, subH, below, pos) + transfer_rightenv!(subenv, below, subH, above, pos) end return envs end function transfer_leftenv!(envs::MultipleEnvironments{<:InfiniteEnvironments}, - above, operator, below, pos::Int) + below, operator, above, pos::Int) for (subH, subenv) in zip(operator, envs.envs) - transfer_leftenv!(subenv, above, subH, below, pos) + transfer_leftenv!(subenv, below, subH, above, pos) end return envs end diff --git a/src/operators/multipliedoperator.jl b/src/operators/multipliedoperator.jl index 77b55fbaa..9b1848196 100644 --- a/src/operators/multipliedoperator.jl +++ b/src/operators/multipliedoperator.jl @@ -52,6 +52,6 @@ Base.:*(op::UntimedOperator, g::Function) = MultipliedOperator(op.op, t -> g(t) function environments(st, x::MultipliedOperator, args...; kwargs...) return environments(st, x.op, args...; kwargs...) end -function recalculate!(envs, above, x::MultipliedOperator, below=above; kwargs...) - return recalculate!(envs, above, x.op, below; kwargs...) +function recalculate!(envs, below, x::MultipliedOperator, above=below; kwargs...) + return recalculate!(envs, below, x.op, above; kwargs...) end diff --git a/test/algorithms.jl b/test/algorithms.jl index dcf05e883..f45319133 100644 --- a/test/algorithms.jl +++ b/test/algorithms.jl @@ -687,6 +687,8 @@ end verbosity = verbosity_conv @testset "mpo * infinite ≈ infinite" begin ψ = InfiniteMPS([ℙ^2, ℙ^2], [ℙ^10, ℙ^10]) + ψ0 = InfiniteMPS([ℙ^2, ℙ^2], [ℙ^12, ℙ^12]) + H = force_planar(repeat(transverse_field_ising(; g=4), 2)) dt = 1e-3 @@ -695,12 +697,12 @@ end W1 = MPSKit.DenseMPO(sW1) W2 = MPSKit.DenseMPO(sW2) - ψ1, _ = approximate(ψ, (sW1, ψ), VOMPS(; verbosity)) - ψ2, _ = approximate(ψ, (W2, ψ), VOMPS(; verbosity)) - ψ3, _ = approximate(ψ, (W1, ψ), IDMRG(; verbosity)) - ψ4, _ = approximate(ψ, (sW2, ψ), IDMRG2(; trscheme=truncdim(20), verbosity)) + ψ1, _ = approximate(ψ0, (sW1, ψ), VOMPS(; verbosity)) + ψ2, _ = approximate(ψ0, (W2, ψ), VOMPS(; verbosity)) + ψ3, _ = approximate(ψ0, (W1, ψ), IDMRG(; verbosity)) + ψ4, _ = approximate(ψ0, (sW2, ψ), IDMRG2(; trscheme=truncdim(12), verbosity)) ψ5, _ = timestep(ψ, H, 0.0, dt, TDVP()) - ψ6 = changebonds(W1 * ψ, SvdCut(; trscheme=truncdim(10))) + ψ6 = changebonds(W1 * ψ, SvdCut(; trscheme=truncdim(12))) @test abs(dot(ψ1, ψ5)) ≈ 1.0 atol = dt @test abs(dot(ψ3, ψ5)) ≈ 1.0 atol = dt