From d35574b6fc3c9a72170440552b8ad4adb6322496 Mon Sep 17 00:00:00 2001 From: LHerviou <98113223+LHerviou@users.noreply.github.com> Date: Wed, 13 Sep 2023 23:45:02 +0200 Subject: [PATCH] Improve `InfiniteMPOMatrix`, rename to `InfiniteBlockMPO` (#77) --- src/ITensorInfiniteMPS.jl | 6 +- src/ITensors.jl | 10 + src/celledvectors.jl | 7 +- src/infiniteblockmpo.jl | 347 +++++++++++++++++++++++++++++++++++ src/infinitempo.jl | 169 +++++++++++++++++ src/infinitempomatrix.jl | 154 ---------------- src/models/fqhe13.jl | 57 ++++-- src/models/models.jl | 127 +++++++++++-- src/subspace_expansion.jl | 2 +- src/vumps_mpo.jl | 125 +++++++------ test/test_iMPOConversions.jl | 117 ++++++++++-- test/test_vumpsmpo.jl | 22 +-- test/test_vumpsmpo_fqhe.jl | 2 +- 13 files changed, 861 insertions(+), 284 deletions(-) create mode 100644 src/infiniteblockmpo.jl delete mode 100644 src/infinitempomatrix.jl diff --git a/src/ITensorInfiniteMPS.jl b/src/ITensorInfiniteMPS.jl index d0ef847d..c1af600a 100644 --- a/src/ITensorInfiniteMPS.jl +++ b/src/ITensorInfiniteMPS.jl @@ -34,9 +34,9 @@ include("itensormap.jl") include("celledvectors.jl") include("abstractinfinitemps.jl") include("infinitemps.jl") -include("infinitempo.jl") include("infinitecanonicalmps.jl") -include("infinitempomatrix.jl") +include("infiniteblockmpo.jl") +include("infinitempo.jl") include("transfermatrix.jl") include("models/models.jl") include("models/fqhe13.jl") @@ -60,7 +60,7 @@ export Cell, InfMPS, InfiniteSum, InfiniteMPO, - InfiniteMPOMatrix, + InfiniteBlockMPO, InfiniteSumLocalOps, ITensorMap, ITensorNetwork, diff --git a/src/ITensors.jl b/src/ITensors.jl index 547df25b..7b60b6f6 100644 --- a/src/ITensors.jl +++ b/src/ITensors.jl @@ -182,3 +182,13 @@ end # TODO: make this definition AbstractMPS # Handle orthogonality center correctly Base.getindex(ψ::MPS, r::UnitRange{Int}) = MPS([ψ[n] for n in r]) + +#TODO Remove if everything is working nicely +#Was still crashing on my laptop after updating ITensors +Base.fill!(::NDTensors.NoData, ::Any) = NDTensors.NoData() + +function ITensors.NDTensors.contraction_output( + A::NDTensors.EmptyTensor, B::NDTensors.DiagBlockSparseTensor, label +) + return NDTensors.EmptyTensor(promote_type(eltype(A), eltype(B)), label) +end diff --git a/src/celledvectors.jl b/src/celledvectors.jl index 16a8ad3d..3686d060 100644 --- a/src/celledvectors.jl +++ b/src/celledvectors.jl @@ -142,6 +142,7 @@ _setindex_cell1!(cv::CelledVector, val, n::Int) = (ITensors.data(cv)[n] = val) function getindex(cv::CelledVector, n::Int) cellₙ = cell(cv, n) siteₙ = cellindex(cv, n) + cellₙ == 1 && return _getindex_cell1(cv, siteₙ) #Avoid unnecessary calls return translatecell(cv.translator, _getindex_cell1(cv, siteₙ), cellₙ - 1) end @@ -169,7 +170,11 @@ getindex(cv::CelledVector, c::Cell) = cv[eachindex(cv, c)] function setindex!(cv::CelledVector, T, n::Int) cellₙ = cell(cv, n) siteₙ = cellindex(cv, n) - _setindex_cell1!(cv, translatecell(cv.translator, T, -(cellₙ - 1)), siteₙ) + if cellₙ == 1 + _setindex_cell1!(cv, T, siteₙ) + else + _setindex_cell1!(cv, translatecell(cv.translator, T, -(cellₙ - 1)), siteₙ) + end return cv end diff --git a/src/infiniteblockmpo.jl b/src/infiniteblockmpo.jl new file mode 100644 index 00000000..5a2bd6c2 --- /dev/null +++ b/src/infiniteblockmpo.jl @@ -0,0 +1,347 @@ + +mutable struct InfiniteBlockMPO <: AbstractInfiniteMPS + data::CelledVector{Matrix{ITensor}} + llim::Int #RealInfinity + rlim::Int #RealInfinity + reverse::Bool +end + +translator(mpo::InfiniteBlockMPO) = mpo.data.translator +data(mpo::InfiniteBlockMPO) = mpo.data + +# TODO better printing? +function Base.show(io::IO, M::InfiniteBlockMPO) + print(io, "$(typeof(M))") + (length(M) > 0) && print(io, "\n") + for i in eachindex(M) + if !isassigned(M, i) + println(io, "#undef") + else + A = M[i] + println(io, "Matrix tensor of size $(size(A))") + for k in 1:size(A, 1), l in 1:size(A, 2) + if !isassigned(A, k + (size(A, 1) - 1) * l) + println(io, "[($k, $l)] #undef") + elseif isempty(A[k, l]) + println(io, "[($k, $l)] empty") + else + println(io, "[($k, $l)] $(inds(A[k, l]))") + end + end + end + end +end + +function getindex(ψ::InfiniteBlockMPO, n::Integer) + return ψ.data[n] +end + +function InfiniteBlockMPO(arrMat::Vector{Matrix{ITensor}}) + return InfiniteBlockMPO(arrMat, 0, size(arrMat, 1), false) +end + +function InfiniteBlockMPO(data::Vector{Matrix{ITensor}}, translator::Function) + return InfiniteBlockMPO(CelledVector(data, translator), 0, size(data, 1), false) +end + +function InfiniteBlockMPO(data::CelledVector{Matrix{ITensor}}, m::Int64, n::Int64) + return InfiniteBlockMPO(data, m, n, false) +end + +function InfiniteBlockMPO(data::CelledVector{Matrix{ITensor}}) + return InfiniteBlockMPO(data, 0, size(data, 1), false) +end + +function ITensors.siteinds(A::InfiniteBlockMPO) + data = [dag(only(filterinds(uniqueinds(A[1][1, 1], A[2][1, 1]); plev=0)))] + for x in 2:(nsites(A) - 1) + append!( + data, + [ + dag( + only(filterinds(uniqueinds(A[x][1, 1], A[x - 1][1, 1], A[x + 1][1, 1]); plev=0)) + ), + ], + ) + end + append!( + data, + [dag(only(filterinds(uniqueinds(A[nsites(A)][1, 1], A[nsites(A) - 1][1, 1]); plev=0)))], + ) + return CelledVector(data, translator(A)) +end + +function ITensors.splitblocks(H::InfiniteBlockMPO) + H = copy(H) + N = nsites(H) + for j in 1:N + for n in 1:length(H) + H[j][n] = splitblocks(H[j][n]) + end + end + return H +end + +function find_all_links(Hm::Matrix{ITensor}) + is = inds(Hm[1, 1]) #site inds + lx, ly = size(Hm) + #We extract the links from the order-3 tensors on the first column and line + #We add dummy indices if there is no relevant indices + ir = only(uniqueinds(Hm[1, 2], is)) + ir0 = Index(ITensors.trivial_space(ir); dir=dir(ir), tags="Link,extra") + il0 = dag(ir0) + left_links = typeof(ir)[] + for x in 1:lx + temp = uniqueinds(Hm[x, 1], is) + if length(temp) == 0 + append!(left_links, [il0]) + elseif length(temp) == 1 + append!(left_links, temp) + else + error("ITensor does not seem to be of the correct order") + end + end + right_links = typeof(ir)[] + for x in 1:lx + temp = uniqueinds(Hm[1, x], is) + if length(temp) == 0 + append!(right_links, [ir0]) + elseif length(temp) == 1 + append!(right_links, temp) + else + error("ITensor does not seem to be of the correct order") + end + end + return left_links, right_links +end + +""" + local_mpo_block_projectors(is::Index; new_tags = tags(is)) + + Build the projectors on the three parts of the itensor used to split a MPO into an InfiniteBlockMPO + More precisely, create projectors on the first dimension, the 2:end-1 and the last dimension of the index + Input: is the Index to split + Output: the triplet of projectors (first, middle, last) + Optional arguments: new_tags: if we want to change the tags of the index. +""" +function local_mpo_block_projectors(is::Index; tags=tags(is)) + old_dim = dim(is) + #Build the local projectors. + #We have to differentiate between the dense and the QN case + #Note that as far as I know, the MPO even dense is always guaranteed to have identities at both corners + #If it is not the case, my construction will not work + top = onehot(dag(is) => 1) + bottom = onehot(dag(is) => old_dim) + if length(is.space) == 1 + new_ind = Index(is.space - 2; tags=tags) + mat = zeros(new_ind.space, is.space) + for x in 1:(new_ind.space) + mat[x, x + 1] = 1 + end + middle = ITensor(copy(mat), new_ind, dag(is)) + else + new_ind = Index(is.space[2:(end - 1)]; dir=dir(is), tags=tags) + middle = ITensors.BlockSparseTensor( + Float64, + undef, + Block{2}[Block(x, x + 1) for x in 1:length(new_ind.space)], + (new_ind, dag(is)), + ) + for x in 1:length(new_ind.space) + dim_block = new_ind.space[x][2] + ITensors.blockview(middle, Block(x, x + 1)) .= diagm(0 => ones(dim_block)) + end + middle = itensor(middle) + end + return top, middle, bottom +end + +""" + local_mpo_blocks(tensor::ITensor, left_ind::Index, right_ind::Index; left_tags = tags(inds[1]), right_tags = tags(inds[2]), ...) + + Converts a 4-legged tensor (coming from an (infinite) MPO) with two site indices and a left and a right leg into a 3 x 3 matrix of ITensor. + We assume the normal form for MPO (before full compression) where the top left and bottom right corners are identity matrices. + The goal is to write the tensor in the form + 1 M_12 M_13 + M_21 M_22 M_23 + M_31 M_32 1 + such that we can then easily compress it. Note that for most of our tensors, the upper triangular part will be 0. + Input: tensor the four leg tensors and the pair of Index (left_ind, right_ind) + Output: the 3x3 matrix of tensors + Optional arguments: left_tags: if we want to change the tags of the left indices. + right_tags: if we want to change the tags of the right indices. +""" +function local_mpo_blocks( + t::ITensor, + left_ind::Index, + right_ind::Index; + left_tags=tags(left_ind), + right_tags=tags(right_ind), + kwargs..., +) + @assert order(t) == 4 + + left_dim = dim(left_ind) + right_dim = dim(right_ind) + #Build the local projectors. + top_left, middle_left, bottom_left = local_mpo_block_projectors(left_ind; tags=left_tags) + top_right, middle_right, bottom_right = local_mpo_block_projectors( + right_ind; tags=right_tags + ) + + matrix = Matrix{ITensor}(undef, 3, 3) + for (idx_left, proj_left) in enumerate([top_left, middle_left, bottom_left]) + for (idx_right, proj_right) in enumerate([top_right, middle_right, bottom_right]) + matrix[idx_left, idx_right] = proj_left * t * proj_right + end + end + return matrix +end + +""" + local_mpo_blocks(t::ITensor, ind::Index; new_tags = tags(ind), position = :first, ...) + + Converts a 3-legged tensor (the extremity of a MPO) with two site indices and one leg into a 3 Vector of ITensor. + We assume the normal form for MPO (before full compression) where the top left and bottom right corners are identity matrices in the bulk. + + Input: tensor the three leg tensors and the index connecting to the rest of the MPO + Output: the 3x1 or 1x3 vector of tensors + Optional arguments: new_tags: if we want to change the tags of the indices. + position: whether we consider the first term in the MPO or the last. +""" +function local_mpo_blocks( + t::ITensor, ind::Index; new_tags=tags(ind), position=:first, kwargs... +) + @assert order(t) == 3 + top, middle, bottom = local_mpo_block_projectors(ind; tags=new_tags) + + if position == :first + vector = Matrix{ITensor}(undef, 1, 3) + else + vector = Matrix{ITensor}(undef, 3, 1) + end + for (idx, proj) in enumerate([top, middle, bottom]) + vector[idx] = proj * t + end + return vector +end + +""" +combineblocks_linkinds_auxiliary(Hcl::InfiniteBlockMPO) + +The workhorse of combineblocks_linkinds. We separated them for ease of maintenance. +Fuse the non-site legs of the infiniteBlockMPO Hcl and the corresponding left L and right R environments. +Preserve the corner structure. +Essentially the inverse of splitblocks. It becomes useful for the very dense MPOs once get after compression sometimes. +Input: Hcl the infiniteBlockMPO +Output: a copy of Hcl fused, and the two array of combiners to apply to left and right environments if needed. +""" +function combineblocks_linkinds_auxiliary(H::InfiniteBlockMPO) + H = copy(H) + N = nsites(H) + for j in 1:(N - 1) + right_dim = size(H[j], 2) + for d in 2:(right_dim - 1) + right_link = only(commoninds(H[j][1, d], H[j + 1][d, 1])) + comb = combiner(right_link; tags=tags(right_link)) + comb_ind = combinedind(comb) + for k in 1:size(H[j], 1) + if isempty(H[j][k, d]) + H[j][k, d] = ITensor(Float64, uniqueinds(H[j][k, d], right_link)..., comb_ind) + else + H[j][k, d] = H[j][k, d] * comb + end + end + for k in 1:size(H[j + 1], 2) + if isempty(H[j + 1][d, k]) + H[j + 1][d, k] = ITensor( + Float64, uniqueinds(H[j + 1][d, k], dag(right_link))..., dag(comb_ind) + ) + else + H[j + 1][d, k] = H[j + 1][d, k] * dag(comb) + end + end + end + end + right_dim = size(H[N], 2) + left_combs = [] + right_combs = [] + for d in 2:(right_dim - 1) + right_link = only(commoninds(H[N][1, d], H[N + 1][d, 1])) + comb = combiner(right_link; tags=tags(right_link)) + comb_ind = combinedind(comb) + comb2 = translatecell(translator(H), comb, -1) + comb_ind2 = translatecell(translator(H), comb_ind, -1) + for k in 1:size(H[N], 1) + if isempty(H[N][k, d]) + H[N][k, d] = ITensor(Float64, uniqueinds(H[N][k, d], right_link)..., comb_ind) + else + H[N][k, d] = H[N][k, d] * comb + end + end + for k in 1:size(H[1], 2) + if isempty(H[1][d, k]) + H[1][d, k] = ITensor( + Float64, + uniqueinds(H[1][d, k], dag(translatecell(translator(H), right_link, -1)))..., + dag(comb_ind2), + ) + else + H[1][d, k] = H[1][d, k] * dag(comb2) + end + end + append!(left_combs, [comb2]) + append!(right_combs, [dag(comb)]) + end + return H, left_combs, right_combs +end + +""" +combineblocks_linkinds(Hcl::InfiniteBlockMPO, left_environment, right_environment) + +Fuse the non-site legs of the infiniteBlockMPO Hcl and the corresponding left L and right R environments. +Preserve the corner structure. +Essentially the inverse of splitblocks. It becomes useful for the very dense MPOs once get after compression sometimes. +Input: Hcl the infiniteBlockMPO, the left environment and right environment +Output: Hcl, left_environment, right_environment the updated MPO and environments +""" +function combineblocks_linkinds(H::InfiniteBlockMPO, left_environment, right_environment) + H, left_combs, right_combs = combineblocks_linkinds_auxiliary(H) + left_environment = copy(left_environment) + for j in 1:length(left_combs) + if isempty(left_environment[j + 1]) + left_environment[j + 1] = ITensor( + uniqueinds(left_environment[j + 1], left_combs[j])..., + uniqueinds(left_combs[j], left_environment[j + 1])..., + ) + else + left_environment[j + 1] = left_environment[j + 1] * left_combs[j] + end + end + right_environment = copy(right_environment) + for j in 1:length(right_combs) + if isempty(right_environment[j + 1]) + right_environment[j + 1] = ITensor( + uniqueinds(right_environment[j + 1], right_combs[j])..., + uniqueinds(right_combs[j], right_environment[j + 1])..., + ) + else + right_environment[j + 1] = right_environment[j + 1] * right_combs[j] + end + end + return H, left_environment, right_environment +end + +""" +combineblocks_linkinds(Hcl::InfiniteBlockMPO) + +Fuse the non-site legs of the infiniteBlockMPO Hcl. +Preserve the corner structure. +Essentially the inverse of splitblocks. It becomes useful for the very dense MPOs once get after compression sometimes. +Input: Hcl the infiniteBlockMPO, +Output: the updated MPO +""" +function combineblocks_linkinds(H::InfiniteBlockMPO) + H, _ = combineblocks_linkinds_auxiliary(H) + return H +end diff --git a/src/infinitempo.jl b/src/infinitempo.jl index db9c44fb..dad99373 100644 --- a/src/infinitempo.jl +++ b/src/infinitempo.jl @@ -13,3 +13,172 @@ end translator(mpo::InfiniteMPO) = mpo.data.translator InfiniteMPO(data::CelledVector{ITensor}) = InfiniteMPO(data, 0, size(data, 1), false) + +# +# Hm should have the form below. Only link indices are shown. +# +# I M-->l=n : : : +# l=n-->M l=n-->M-->l=n : : : +# l=n-1-->M l=n-1-->M-->l=n ... ; : : +# : : : : : +# l=1-->M l=1->M-->l=n ... l=1-->M-->l=2 l=1-->M-->l=1 l=1-->M +# M M-->l=n ... M-->l=2 M-->l=1 I +# +# We no longer assume any specific form of the matrix of tensor, and we squash all elements together +# This is facilutated by making all elements of Hm into order(4) tensors by adding dummy Dw=1 indices. +# We use on line then columns ITensors.directsum() to join all the blocks into the final ITensor. +# This code should be 100% dense/blocks-sparse agnostic. +# +function cat_to_itensor(Hm::Matrix{ITensor}) + lx, ly = size(Hm) + T = eltype(Hm[1, 1]) + left_links, right_links = find_all_links(Hm) + #Start by fusing lines together + Ls = [] + new_rs = [] + for x in 1:lx + #Formatting issues force me to write this as a loop + input = [Hm[x, 1] * onehot(T, right_links[1] => 1) => right_links[1]] + for y in 2:(ly - 1) + append!(input, [Hm[x, y] => right_links[y]]) + end + append!(input, [Hm[x, ly] * onehot(T, right_links[ly] => 1) => right_links[ly]]) + H, ir = directsum(input...; tags="Link, right") + if x == 1 + append!(new_rs, [ir]) + else + replaceinds!(H, [ir], [new_rs[1]]) + end + append!(Ls, [H]) + end + input = [Ls[1] * onehot(T, left_links[1] => 1) => left_links[1]] + for x in 2:(lx - 1) + append!(input, [Ls[x] => left_links[x]]) + end + append!(input, [Ls[lx] * onehot(T, left_links[lx] => 1) => left_links[lx]]) + H, new_l = directsum(input...; tags="Link, left") + return H, new_l, new_rs[1] +end +# +# Hm is the InfiniteBlockMPO +# Hlrs is an array of {ITensor,Index,Index}s, one for each site in the unit cell. +# Hi is a CelledVector of ITensors. +# +function InfiniteMPO(Hm::InfiniteBlockMPO) + Hlrs = cat_to_itensor.(Hm) #return an array of {ITensor,Index,Index} + # + # Unpack the array of tuples into three arrays. And also get an array site indices. + # + Hi = CelledVector([Hlr[1] for Hlr in Hlrs], translator(Hm)) + ils = CelledVector([Hlr[2] for Hlr in Hlrs], translator(Hm)) + irs = CelledVector([Hlr[3] for Hlr in Hlrs], translator(Hm)) + s = siteinds(Hm) + # + # Create new tags with proper cell and link numbers. Also daisy chain + # all the indices so right index at j = dag(left index at j+1) + # + for j in 1:nsites(Hm) + newTag = "Link,c=$(getcell(s[j])),l=$(getsite(s[j]))" + ir = replacetags(irs[j], tags(irs[j]), newTag) #new right index + Hi[j] = replaceinds(Hi[j], [irs[j]], [ir]) + Hi[j + 1] = replaceinds(Hi[j + 1], [ils[j + 1]], [dag(ir)]) + end + return InfiniteMPO(Hi) +end + +""" + combineblocks_linkinds_auxiliary(Hcl::InfiniteMPO) + + The workhorse of combineblocks_linkinds. We separated them for ease of maintenance. + Fuse the non-site legs of the infiniteMPO Hcl and the corresponding left L and right R environments. + Preserve the corner structure. + Essentially the inverse of splitblocks. It becomes useful for the very dense MPOs once get after compression sometimes. + Input: Hcl the infinite MPO + Output: Hcl, the left and right combiners to use on environments +""" +function combineblocks_linkinds_auxiliary(Hcl::InfiniteMPO) + Hcl = copy(Hcl) + N = nsites(Hcl) + for j in 1:(N - 1) + right_link = only(commoninds(Hcl[j], Hcl[j + 1])) + top, middle, bottom = local_mpo_block_projectors(right_link) + if length(right_link.space) == 1 + top = top * onehot(Index(1) => 1) + bottom = bottom * onehot(Index(1) => 1) + else + top = top * onehot(Index(QN() => 1) => 1) + bottom = bottom * onehot(Index(QN() => 1) => 1) + end + cb = combiner( + only(uniqueinds(middle, dag(right_link))); dir=dir(right_link), tags=tags(right_link) + ) + middle = cb * middle + comb, comb_ind = directsum( + top => uniqueinds(top, dag(right_link)), + middle => uniqueinds(middle, dag(right_link)), + bottom => uniqueinds(bottom, dag(right_link)); + tags=[tags(right_link)], + ) + Hcl[j] = Hcl[j] * comb + Hcl[j + 1] = dag(comb) * Hcl[j + 1] + end + right_link = only(commoninds(Hcl[N], Hcl[N + 1])) + right_link2 = translatecell(translator(Hcl), right_link, -1) + top, middle, bottom = local_mpo_block_projectors(right_link) + if length(right_link.space) == 1 + top = top * onehot(Index(1) => 1) + bottom = bottom * onehot(Index(1) => 1) + else + top = top * onehot(Index(QN() => 1) => 1) + bottom = bottom * onehot(Index(QN() => 1) => 1) + end + cb = combiner( + only(uniqueinds(middle, dag(right_link))); dir=dir(right_link), tags=tags(right_link) + ) + middle = cb * middle + comb, comb_ind = directsum( + top => uniqueinds(top, dag(right_link)), + middle => uniqueinds(middle, dag(right_link)), + bottom => uniqueinds(bottom, dag(right_link)); + tags=[tags(right_link)], + ) + comb2 = translatecell(translator(Hcl), comb, -1) + comb2_ind = translatecell(translator(Hcl), comb2, -1) + Hcl[N] = Hcl[N] * comb + Hcl[1] = dag(comb2) * Hcl[1] + return Hcl, comb2, dag(comb) +end + +""" + combineblocks_linkinds(Hcl::InfiniteMPO, left_environment::ITensor, right_environment::ITensor) + + Fuse the non-site legs of the infiniteMPO Hcl and the corresponding left and right environments. + Preserve the corner structure. + Essentially the inverse of splitblocks. It becomes useful for the very dense MPOs once get after compression sometimes. + Hcl is modified on site, L and R are not. + Input: Hcl the infinite MPO + left_environment the left environment (an ITensor) + right_environment the right environment (an ITensor) + Output: Hcl, newL, newR the updated MPO and environments +""" +function combineblocks_linkinds( + Hcl::InfiniteMPO, left_environment::ITensor, right_environment::ITensor +) + Hcl, left_comb, right_comb = combineblocks_linkinds_auxiliary(Hcl) + return Hcl, left_environment * left_comb, right_comb * right_environment +end + +""" + combineblocks_linkinds(Hcl::InfiniteMPO) + + Fuse the non-site legs of the infiniteMPO Hcl. + Preserve the corner structure. + Essentially the inverse of splitblocks. It becomes useful for the very dense MPOs once get after compression sometimes. + Hcl is modified on site, L and R are not. + Input: Hcl the infinite MPO + Output: Hcl, +""" +function combineblocks_linkinds(Hcl::InfiniteMPO) + Hcl, = combineblocks_linkinds_auxiliary(Hcl) + return Hcl +end diff --git a/src/infinitempomatrix.jl b/src/infinitempomatrix.jl deleted file mode 100644 index 29ec9667..00000000 --- a/src/infinitempomatrix.jl +++ /dev/null @@ -1,154 +0,0 @@ - -mutable struct InfiniteMPOMatrix <: AbstractInfiniteMPS - data::CelledVector{Matrix{ITensor}} - llim::Int #RealInfinity - rlim::Int #RealInfinity - reverse::Bool -end - -translator(mpo::InfiniteMPOMatrix) = mpo.data.translator -data(mpo::InfiniteMPOMatrix) = mpo.data - -# TODO better printing? -function Base.show(io::IO, M::InfiniteMPOMatrix) - print(io, "$(typeof(M))") - (length(M) > 0) && print(io, "\n") - for i in eachindex(M) - if !isassigned(M, i) - println(io, "#undef") - else - A = M[i] - println(io, "Matrix tensor of size $(size(A))") - for k in 1:size(A)[1], l in 1:size(A)[2] - if !isassigned(A, k + (size(A)[1] - 1) * l) - println(io, "[($k, $l)] #undef") - elseif isempty(A[k, l]) - println(io, "[($k, $l)] empty") - else - println(io, "[($k, $l)] $(inds(A[k, l]))") - end - end - end - end -end - -function getindex(ψ::InfiniteMPOMatrix, n::Integer) - return ψ.data[n] -end - -function InfiniteMPOMatrix(arrMat::Vector{Matrix{ITensor}}) - return InfiniteMPOMatrix(arrMat, 0, size(arrMat)[1], false) -end - -function InfiniteMPOMatrix(data::Vector{Matrix{ITensor}}, translator::Function) - return InfiniteMPOMatrix(CelledVector(data, translator), 0, size(data)[1], false) -end - -# -# Hm should have the form below. Only link indices are shown. -# -# I 0 0 0 0 -# M-->l=n 0 0 0 0 -# 0 l=n-->M-->l=n-1 ... 0 0 0 -# : : : : : -# 0 0 ... l=2-->M-->l=1 0 0 -# 0 0 ... 0 l=1-->M I -# -# We need to capture the corner I's and the sub-diagonal Ms, and paste them together into an ITensor. -# This is facilutated by making all elements of Hm into order(4) tensors by adding dummy Dw=1 indices. -# We ITensors.directsum() to join all the blocks into the final ITensor. -# This code should be 100% dense/blocks-sparse agnostic. -# -function cat_to_itensor(Hm::Matrix{ITensor}, inds_array)::Tuple{ITensor,Index,Index} - lx, ly = size(Hm) - @assert lx == ly - # - # Extract the sub diagonal - # - Ms = map(n -> Hm[lx - n + 1, ly - n], 1:(lx - 1)) #Get the sub diagonal into an array. - N = length(Ms) - @assert N == lx - 1 - # - # Convert edge Ms to order 4 ITensors using the dummy index. - # - il0 = inds_array[1][1] - T = eltype(Ms[1]) - Ms[1] *= onehot(T, il0 => 1) - Ms[N] *= onehot(T, dag(il0) => 1) #We don't need distinct index IDs for this. directsum makes all new index anyway. - # - # Start direct sums to build the single ITensor. Order is critical to get the desired form. - # - H = Hm[1, 1] * onehot(T, il0' => 1) * onehot(T, dag(il0) => 1) #Bootstrap with the I op in the top left of Hm - H, il = directsum(H => il0', Ms[N] => inds_array[N][1]) #1-D directsum to put M[N] directly below the I op - ir = dag(il0) #setup recusion. - for i in (N - 1):-1:1 #2-D directsums to place new blocks below and to the right. - H, (il, ir) = directsum( - H => (il, ir), Ms[i] => inds_array[i]; tags=["Link,left", "Link,right"] - ) - end - IN = Hm[N + 1, N + 1] * onehot(T, dag(il0) => 1) * onehot(T, il => dim(il)) #I op in the bottom left of Hm - H, ir = directsum(H => ir, IN => il0; tags="Link,right") #1-D directsum to the put I in the bottom row, to the right of M[1] - - return H, il, ir -end - -function find_all_links(Hms::InfiniteMPOMatrix) - is = inds(Hms[1][1, 1]) #site inds - ir = noncommonind(Hms[1][2, 1], is) #This op should have one link index pointing to the right. - # - # Set up return array of 2-tuples - # - indexT = typeof(ir) - TupleT = NTuple{2,indexT} - inds_array = Vector{TupleT}[] - # - # Make a dummy index - # - il0 = Index(ITensors.trivial_space(ir); dir=dir(dag(ir)), tags="Link,l=0") - # - # Loop over sites and nonzero matrix elements which are linked into the next - # and previous iMPOMatrix. - # - for n in 1:nsites(Hms) - Hm = Hms[n] - inds_n = TupleT[] - lx, ly = size(Hm) - @assert lx == ly - for iy in (ly - 1):-1:1 - ix = iy + 1 - il = ix < lx ? commonind(Hm[ix, iy], Hms[n - 1][ix + 1, iy + 1]) : dag(il0) - ir = iy > 1 ? commonind(Hm[ix, iy], Hms[n + 1][ix - 1, iy - 1]) : il0 - push!(inds_n, (il, ir)) - end - push!(inds_array, inds_n) - end - return inds_array -end - -# -# Hm is the InfiniteMPOMatrix -# Hlrs is an array of {ITensor,Index,Index}s, one for each site in the unit cell. -# Hi is a CelledVector of ITensors. -# -function InfiniteMPO(Hm::InfiniteMPOMatrix) - inds_array = find_all_links(Hm) - Hlrs = cat_to_itensor.(Hm, inds_array) #return an array of {ITensor,Index,Index} - # - # Unpack the array of tuples into three arrays. And also get an array site indices. - # - Hi = CelledVector([Hlr[1] for Hlr in Hlrs], translator(Hm)) - ils = CelledVector([Hlr[2] for Hlr in Hlrs], translator(Hm)) - irs = CelledVector([Hlr[3] for Hlr in Hlrs], translator(Hm)) - site_inds = [commoninds(Hm[j][1, 1], Hm[j][end, end])[1] for j in 1:nsites(Hm)] - # - # Create new tags with proper cell and link numbers. Also daisy chain - # all the indices so right index at j = dag(left index at j+1) - # - for j in 1:nsites(Hm) - newTag = "Link,c=$(getcell(site_inds[j])),l=$(getsite(site_inds[j]))" - ir = replacetags(irs[j], tags(irs[j]), newTag) #new right index - Hi[j] = replaceinds(Hi[j], [irs[j]], [ir]) - Hi[j + 1] = replaceinds(Hi[j + 1], [ils[j + 1]], [dag(ir)]) - end - return InfiniteMPO(Hi) -end diff --git a/src/models/fqhe13.jl b/src/models/fqhe13.jl index 40d87e6e..79a7fd9d 100644 --- a/src/models/fqhe13.jl +++ b/src/models/fqhe13.jl @@ -1,12 +1,43 @@ +function ITensors.space(::SiteType"FermionK", pos::Int; p=1, q=1, conserve_momentum=true) + if !conserve_momentum + return [QN("Nf", -p) => 1, QN("Nf", q - p) => 1] + else + return [ + QN(("Nf", -p), ("NfMom", -p * pos)) => 1, + QN(("Nf", q - p), ("NfMom", (q - p) * pos)) => 1, + ] + end +end + +# Forward all op definitions to Fermion +function ITensors.op!(Op::ITensor, opname::OpName, ::SiteType"FermionK", s::Index...) + return ITensors.op!(Op, opname, SiteType("Fermion"), s...) +end + +function ITensors.has_fermion_string(opname::OpName, ::SiteType"FermionK") + return has_fermion_string(opname, SiteType("Fermion")) +end + function unit_cell_terms( - ::Model"fqhe_2b_pot"; Ly::Float64, Vs::Array{Float64,1}, prec::Float64 + ::Model"fqhe_2b_pot"; Ly::Float64, Vs::Array{Float64,1}, prec::Float64, chemical_shift=0.0 ) - rough_N = round(Int64, 2 * Ly) - coeff = build_two_body_coefficient_pseudopotential(; N_phi=rough_N, Ly=Ly, Vs=Vs) - opt = optimize_coefficients(coeff; prec=prec) - opt = filter_optimized_Hamiltonian_by_first_site(opt) - #sorted_opt = sort_by_configuration(opt); - return generate_Hamiltonian(opt) + rough_N = round(Int64, 2 * Ly) - 2 + test = round(Int64, 2 * Ly) - 2 + while rough_N <= test + rough_N = test + 2 + coeff = build_two_body_coefficient_pseudopotential(; N_phi=rough_N, Ly=Ly, Vs=Vs) + opt = optimize_coefficients(coeff; prec=prec) + opt = filter_optimized_Hamiltonian_by_first_site(opt) + if haskey(opt, ["N", 1]) + opt[["N", 1]] += chemical_shift + elseif chemical_shift != 0 + opt[["N", 1]] = chemical_shift + end + test = check_max_range_optimized_Hamiltonian(opt) + if rough_N > test + return generate_Hamiltonian(opt) + end + end end #Please contact Loic Herviou before using this part of the code for production @@ -270,16 +301,10 @@ function optimize_coefficients(coeff::Dict; prec=1e-12) if abs(v) < prec continue end - if length(ke) == 4 - name = ["Cdag", "Cdag", "C", "C"] - elseif length(ke) == 6 - name = ["Cdag", "Cdag", "Cdag", "C", "C", "C"] - elseif length(ke) == 8 - name = ["Cdag", "Cdag", "Cdag", "Cdag", "C", "C", "C", "C"] - else - println("Optimization not implemented") - continue + if mod(length(ke), 2) == 1 + error("Odd number of operators is not implemented") end + name = vcat(fill("Cdag", length(ke) ÷ 2), fill("C", length(ke) ÷ 2)) k = Base.copy(ke) sg = get_perm!(k, name) filter_op!(k, name) diff --git a/src/models/models.jl b/src/models/models.jl index 6646d03e..206defd4 100644 --- a/src/models/models.jl +++ b/src/models/models.jl @@ -52,7 +52,7 @@ function shift_sites(term::Scaled{C,Prod{Op}}, shift::Int) where {C} end # Shift the sites of the terms of the OpSum by shift. -# By default, it shifts +# By default, it shifts function shift_sites(opsum::OpSum, shift::Int) shifted_opsum = OpSum() for o in ITensors.terms(opsum) @@ -87,21 +87,27 @@ function InfiniteMPO(model::Model, s::CelledVector; kwargs...) end function InfiniteMPO(model::Model, s::CelledVector, translator::Function; kwargs...) - return InfiniteMPO(InfiniteMPOMatrix(model, s, translator; kwargs...)) + return InfiniteMPO(InfiniteBlockMPO(model, s, translator; kwargs...)) end -function InfiniteMPOMatrix(model::Model, s::CelledVector; kwargs...) - return InfiniteMPOMatrix(model, s, translator(s); kwargs...) +function InfiniteBlockMPO(model::Model, s::CelledVector; kwargs...) + return InfiniteBlockMPO(model, s, translator(s); kwargs...) end -function InfiniteMPOMatrix(model::Model, s::CelledVector, translator::Function; kwargs...) +function InfiniteBlockMPO(model::Model, s::CelledVector, translator::Function; kwargs...) N = length(s) temp_H = InfiniteSum{MPO}(model, s; kwargs...) - range_H = nrange(temp_H)[1] - ls = CelledVector([Index(1, "Link,c=1,n=$n") for n in 1:N], translator) + range_H = maximum(nrange(temp_H)) #Should be improved + ls = CelledVector( + [Index(ITensors.trivial_space(s[n]), "Link,c=1,n=$n") for n in 1:N], translator + ) mpos = [Matrix{ITensor}(undef, 1, 1) for i in 1:N] for j in 1:N - Hmat = fill(op("Zero", s[j]), range_H + 1, range_H + 1) + #For type stability + + Hmat = fill( + ITensor(eltype(temp_H[1][1]), dag(s[j]), prime(s[j])), range_H + 1, range_H + 1 + ) identity = op("Id", s[j]) Hmat[1, 1] = identity Hmat[end, end] = identity @@ -110,14 +116,110 @@ function InfiniteMPOMatrix(model::Model, s::CelledVector, translator::Function; if isnothing(idx) Hmat[range_H + 1 - n, range_H - n] = identity else - Hmat[range_H + 1 - n, range_H - n] = temp_H[j - n][idx]#replacetags(linkinds, temp_H[j - n][idx], "Link, l=$n", tags(ls[j-1])) + #Here, we split the local tensor into its different blocks + T = eltype(temp_H[j - n][idx]) + if n == 0 + ind_to_change = only(commoninds(temp_H[j - n][idx], temp_H[j - n][idx + 1])) + temp_mat = local_mpo_blocks( + temp_H[j - n][idx], ind_to_change; position=:first, new_tags=tags(ls[j]) + ) + elseif n == range_H - 1 + ind_to_change = only(commoninds(temp_H[j - n][idx], temp_H[j - n][idx - 1])) + temp_mat = local_mpo_blocks( + temp_H[j - n][idx], ind_to_change; position=:last, new_tags=tags(ls[j - 1]) + ) + else + left_dir = only(commoninds(temp_H[j - n][idx], temp_H[j - n][idx - 1])) + right_dir = only(commoninds(temp_H[j - n][idx], temp_H[j - n][idx + 1])) + temp_mat = local_mpo_blocks( + temp_H[j - n][idx], + left_dir, + right_dir; + left_tags=tags(ls[j - 1]), + right_tags=tags(ls[j]), + ) + end + if size(temp_mat) == (3, 3) + @assert iszero(temp_mat[1, 2]) + @assert iszero(temp_mat[1, 3]) + @assert iszero(temp_mat[2, 3]) + @assert temp_mat[1, 1] == identity + @assert temp_mat[3, 3] == identity + Hmat[range_H + 1 - n, range_H - n] = temp_mat[2, 2] + Hmat[end, range_H - n] = temp_mat[3, 2] + Hmat[range_H + 1 - n, 1] = temp_mat[2, 1] + elseif size(temp_mat) == (1, 3) + @assert n == 0 + @assert temp_mat[1, 3] == identity + #@assert isempty(temp_mat[1, 1]) || iszero(temp_mat[1, 1]) + Hmat[range_H + 1 - n, range_H - n] = temp_mat[1, 2] + Hmat[range_H + 1 - n, 1] = temp_mat[1, 1] + elseif size(temp_mat) == (3, 1) + @assert (range_H - n) == 1 + @assert temp_mat[1, 1] == identity + #@assert isempty(temp_mat[3, 1]) || iszero(temp_mat[3, 1]) + Hmat[range_H + 1 - n, range_H - n] = temp_mat[2, 1] + Hmat[end, range_H - n] += temp_mat[3, 1] #LH This should do nothing #TODO check + else + error("Unexpected matrix form") + end end end mpos[j] = Hmat #mpos[j] += dense(Hmat) * setelt(ls[j-1] => total_dim) * setelt(ls[j] => total_dim) end - #return mpos - return InfiniteMPOMatrix(mpos, translator) + #unify_indices and add virtual indices to the empty tensors + mpos = InfiniteBlockMPO(mpos, translator) + for x in 1:N + sp = prime(s[x]) + sd = dag(s[x]) + left_inds = [ + only(uniqueinds(mpos[x][j, 1], mpos[x][1, 1])) for j in 2:(size(mpos[x], 1) - 1) + ] + right_inds = [ + only(uniqueinds(mpos[x][end, j], mpos[x][1, 1])) for j in 2:(size(mpos[x], 2) - 1) + ] + if x == N + new_right_inds = [ + dag(only(uniqueinds(mpos[1][j, 1], mpos[1][1, 1]))) for + j in 2:(size(mpos[1], 1) - 1) + ] + for j in 1:length(new_right_inds) + new_right_inds[j] = translatecell(translator, new_right_inds[j], 1) + end + else + new_right_inds = [ + dag(only(uniqueinds(mpos[x + 1][j, 1], mpos[x + 1][1, 1]))) for + j in 2:(size(mpos[x], 2) - 1) + ] + end + for j in 2:(size(mpos[x], 1) - 1) + for k in 2:(size(mpos[x], 2) - 1) + if isempty(mpos[x][j, k]) + mpos[x][j, k] = ITensor(left_inds[j - 1], sd, sp, new_right_inds[k - 1]) + else + replaceinds!(mpos[x][j, k], right_inds[k - 1] => new_right_inds[k - 1]) + end + end + end + for j in [1, size(mpos[x], 1)] + for k in 2:(size(mpos[x], 2) - 1) + if isempty(mpos[x][j, k]) + mpos[x][j, k] = ITensor(sd, sp, new_right_inds[k - 1]) + else + replaceinds!(mpos[x][j, k], right_inds[k - 1] => new_right_inds[k - 1]) + end + end + end + for j in 2:(size(mpos[x], 1) - 1) + for k in [1, size(mpos[x], 2)] + if isempty(mpos[x][j, k]) + mpos[x][j, k] = ITensor(left_inds[j - 1], sd, sp) + end + end + end + end + return mpos end function ITensors.MPO(model::Model, s::Vector{<:Index}; kwargs...) @@ -125,7 +227,7 @@ function ITensors.MPO(model::Model, s::Vector{<:Index}; kwargs...) return splitblocks(linkinds, MPO(opsum, s)) end -translatecell(translator, opsum::OpSum, n::Integer) = translator(opsum, n) +translatecell(translator::Function, opsum::OpSum, n::Integer) = translator(opsum, n) function infinite_terms(model::Model; kwargs...) # An `OpSum` storing all of the terms in the @@ -230,7 +332,6 @@ function ITensors.ITensor(model::Model, s::CelledVector, n::Int64; kwargs...) opsum = shift_sites(opsum, -first_site(opsum) + 1) site_range = n:(n + last_site(opsum) - 1) return contract(MPO(opsum, [s[j] for j in site_range])) - # Deprecated version # return contract(MPO(model, s, n; kwargs...)) end diff --git a/src/subspace_expansion.jl b/src/subspace_expansion.jl index 5d2eb202..f2c9a231 100644 --- a/src/subspace_expansion.jl +++ b/src/subspace_expansion.jl @@ -118,7 +118,7 @@ function generate_twobody_nullspace( end function generate_twobody_nullspace( - ψ::InfiniteCanonicalMPS, H::InfiniteMPOMatrix, b::Tuple{Int,Int}; atol=1e-2 + ψ::InfiniteCanonicalMPS, H::InfiniteBlockMPO, b::Tuple{Int,Int}; atol=1e-2 ) n_1, n_2 = b L, _ = left_environment(H, ψ) diff --git a/src/vumps_mpo.jl b/src/vumps_mpo.jl index e632278e..4a7c5613 100644 --- a/src/vumps_mpo.jl +++ b/src/vumps_mpo.jl @@ -1,7 +1,7 @@ # Struct for use in linear system solver struct AOᴸ ψ::InfiniteCanonicalMPS - H::InfiniteMPOMatrix + H::InfiniteBlockMPO n::Int end @@ -28,23 +28,37 @@ function (A::AOᴸ)(x) return xT - xR end +function initialize_left_environment( + H::InfiniteBlockMPO, ψ::InfiniteCanonicalMPS, n::Int64; init_last=false +) + dₕ = size(H[n + 1], 1) + sit = inds(H[n + 1][1, 1]) + link = commonind(ψ.AL[n], ψ.AL[n + 1]) + Ls = Vector{ITensor}(undef, dₕ) + Ls[1] = ITensor(Float64, link, dag(prime(link))) + if init_last + Ls[end] = denseblocks(δ(link, dag(prime(link)))) + else + Ls[end] = ITensor(Float64, link, dag(prime(link))) + end + for j in 2:(dₕ - 1) + mpo_link = only(uniqueinds(H[n + 1][j, 1], sit)) + Ls[j] = ITensor(Float64, dag(mpo_link), link, dag(prime(link))) + end + return Ls +end + function apply_local_left_transfer_matrix( - Lstart::Vector{ITensor}, H::InfiniteMPOMatrix, ψ::InfiniteCanonicalMPS, n_1::Int64 + Lstart::Vector{ITensor}, H::InfiniteBlockMPO, ψ::InfiniteCanonicalMPS, n_1::Int64 ) dₕ = length(Lstart) ψ′ = dag(ψ)' - Ltarget = Vector{ITensor}(undef, size(H[n_1])[1]) + Ltarget = initialize_left_environment(H, ψ, n_1; init_last=false) for j in 1:dₕ - init = false for k in reverse(j:dₕ) - if !isempty(H[n_1][k, j]) && isassigned(Lstart, k) && !isempty(Lstart[k]) - if isassigned(Ltarget, j) && init - Ltarget[j] += Lstart[k] * ψ.AL[n_1] * ψ′.AL[n_1] * H[n_1][k, j] - else - Ltarget[j] = Lstart[k] * ψ.AL[n_1] * ψ′.AL[n_1] * H[n_1][k, j] - init = true - end + if !isempty(H[n_1][k, j]) && !isempty(Lstart[k]) + Ltarget[j] .+= Lstart[k] * ψ.AL[n_1] * H[n_1][k, j] * ψ′.AL[n_1] end end end @@ -53,25 +67,18 @@ end # apply the left transfer matrix at position n1 to the vector Lstart considering it at position m, adding to Ltarget function apply_local_left_transfer_matrix( - Lstart::ITensor, - m::Int64, - H::InfiniteMPOMatrix, - ψ::InfiniteCanonicalMPS, - n_1::Int64; - reset=true, + Lstart::ITensor, m::Int64, H::InfiniteBlockMPO, ψ::InfiniteCanonicalMPS, n_1::Int64; ) - Ltarget = Vector{ITensor}(undef, size(H[n_1])[1]) + Ltarget = initialize_left_environment(H, ψ, n_1; init_last=false) for j in 1:m - if !isempty(H[n_1][m, j]) - Ltarget[j] = Lstart * ψ.AL[n_1] * H[n_1][m, j] * dag(prime(ψ.AL[n_1])) - end + Ltarget[j] = Lstart * ψ.AL[n_1] * H[n_1][m, j] * dag(prime(ψ.AL[n_1])) end return Ltarget end #apply the left transfer matrix n1:n1+nsites(ψ)-1 function apply_left_transfer_matrix( - Lstart::ITensor, m::Int64, H::InfiniteMPOMatrix, ψ::InfiniteCanonicalMPS, n_1::Int64 + Lstart::ITensor, m::Int64, H::InfiniteBlockMPO, ψ::InfiniteCanonicalMPS, n_1::Int64 ) Ltarget = apply_local_left_transfer_matrix(Lstart, m, H, ψ, n_1) for j in 1:(nsites(ψ) - 1) @@ -82,7 +89,7 @@ end # Also input C bond matrices to help compute the right fixed points # of ψ (R ≈ C * dag(C)) -function left_environment(H::InfiniteMPOMatrix, ψ::InfiniteCanonicalMPS; tol=1e-10) +function left_environment(H::InfiniteBlockMPO, ψ::InfiniteCanonicalMPS; tol=1e-10) N = nsites(H) @assert N == nsites(ψ) @@ -100,10 +107,10 @@ function left_environment(H::InfiniteMPOMatrix, ψ::InfiniteCanonicalMPS; tol=1e eₗ = [0.0] dₕ = size(H[1])[1] - Ls = [Vector{ITensor}(undef, dₕ) for j in 1:N] + #Ls = [Vector{ITensor}(undef, dₕ) for j in 1:N] + Ls = [initialize_left_environment(H, ψ, j; init_last=true) for j in 1:N] #Building the L vector for n_1 = 1 # TM is 2 3 ... N 1 - Ls[1][end] = δˡ(1) # exact by construction localR = ψ.C[1] * δʳ(1) * ψ′.C[1] #to revise for n in reverse(1:(dₕ - 1)) temp_Ls = apply_left_transfer_matrix( @@ -144,7 +151,7 @@ end # Struct for use in linear system solver struct AOᴿ ψ::InfiniteCanonicalMPS - H::InfiniteMPOMatrix + H::InfiniteBlockMPO n::Int end @@ -171,43 +178,36 @@ function (A::AOᴿ)(x) return xT - xR end -# apply the left transfer matrix at position n1 to the vector Lstart, replacing Ltarget -function apply_local_right_transfer_matrix!( - Lstart::Vector{ITensor}, H::InfiniteMPOMatrix, ψ::InfiniteCanonicalMPS, n_1::Int64 +function initialize_right_environment( + H::InfiniteBlockMPO, ψ::InfiniteCanonicalMPS, n::Int64; init_first=false ) - dₕ = length(Lstart) - ψ′ = dag(ψ)' - for j in reverse(1:dₕ) - init = false - for k in reverse(1:j) - if !isempty(H[n_1][j, k]) && isassigned(Lstart, k) && !isempty(Lstart[k]) - if isassigned(Lstart, j) && init - Lstart[j] += Lstart[k] * ψ.AR[n_1] * H[n_1][j, k] * ψ′.AR[n_1] - else - Lstart[j] = Lstart[k] * ψ.AR[n_1] * H[n_1][j, k] * ψ′.AR[n_1] - init = true - end - end - end + dₕ = size(H[n - 1], 2) + sit = inds(H[n - 1][1, 1]) + link = commonind(ψ.AR[n], ψ.AR[n - 1]) + Rs = Vector{ITensor}(undef, dₕ) + Rs[end] = ITensor(Float64, link, dag(prime(link))) + if init_first + Rs[1] = denseblocks(δ(link, dag(prime(link)))) + else + Rs[1] = ITensor(Float64, link, dag(prime(link))) + end + for j in 2:(dₕ - 1) + mpo_link = only(uniqueinds(H[n - 1][1, j], sit)) + Rs[j] = ITensor(Float64, dag(mpo_link), link, dag(prime(link))) end + return Rs end function apply_local_right_transfer_matrix( - Lstart::Vector{ITensor}, H::InfiniteMPOMatrix, ψ::InfiniteCanonicalMPS, n_1::Int64 + Lstart::Vector{ITensor}, H::InfiniteBlockMPO, ψ::InfiniteCanonicalMPS, n_1::Int64 ) dₕ = length(Lstart) ψ′ = dag(ψ)' - Ltarget = Vector{ITensor}(undef, size(H[n_1])[1]) + Ltarget = initialize_right_environment(H, ψ, n_1) for j in reverse(1:dₕ) - init = false for k in reverse(1:j) if !isempty(H[n_1][j, k]) && isassigned(Lstart, k) && !isempty(Lstart[k]) - if isassigned(Ltarget, j) && init - Ltarget[j] += Lstart[k] * ψ.AR[n_1] * H[n_1][j, k] * ψ′.AR[n_1] - else - Ltarget[j] = Lstart[k] * ψ.AR[n_1] * H[n_1][j, k] * ψ′.AR[n_1] - init = true - end + Ltarget[j] += Lstart[k] * ψ.AR[n_1] * H[n_1][j, k] * ψ′.AR[n_1] end end end @@ -218,17 +218,17 @@ end function apply_local_right_transfer_matrix( Lstart::ITensor, m::Int64, - H::InfiniteMPOMatrix, + H::InfiniteBlockMPO, ψ::InfiniteCanonicalMPS, n_1::Int64; reset=true, ) dₕ = size(H[n_1])[1] ψ′ = dag(prime(ψ.AR[n_1])) - Ltarget = Vector{ITensor}(undef, dₕ) + Ltarget = initialize_right_environment(H, ψ, n_1) for j in m:dₕ if !isempty(H[n_1][j, m]) - Ltarget[j] = Lstart * ψ.AR[n_1] * H[n_1][j, m] * ψ′ #TODO optimize + Ltarget[j] = Lstart * ψ.AR[n_1] * H[n_1][j, m] * ψ′ end end return Ltarget @@ -236,7 +236,7 @@ end #apply the right transfer matrix n1:n1+nsites(ψ)-1 function apply_right_transfer_matrix( - Lstart::ITensor, m::Int64, H::InfiniteMPOMatrix, ψ::InfiniteCanonicalMPS, n_1::Int64 + Lstart::ITensor, m::Int64, H::InfiniteBlockMPO, ψ::InfiniteCanonicalMPS, n_1::Int64 ) Ltarget = apply_local_right_transfer_matrix(Lstart, m, H, ψ, n_1) for j in 1:(nsites(ψ) - 1) @@ -245,7 +245,7 @@ function apply_right_transfer_matrix( return Ltarget end -function right_environment(H::InfiniteMPOMatrix, ψ::InfiniteCanonicalMPS; tol=1e-10) +function right_environment(H::InfiniteBlockMPO, ψ::InfiniteCanonicalMPS; tol=1e-10) N = nsites(H) @assert N == nsites(ψ) @@ -258,10 +258,9 @@ function right_environment(H::InfiniteMPOMatrix, ψ::InfiniteCanonicalMPS; tol=1 eᵣ = [0.0] dₕ = size(H[1])[1] - Rs = [Vector{ITensor}(undef, dₕ) for j in 1:N] - #Building the L vector for n_1 = 1 + Rs = [initialize_right_environment(H, ψ, j; init_first=true) for j in 1:N] + #Building the R vector for n_1 = 1 # TM is 2-N 3-N ... 0 - Rs[1][1] = δʳ(0) localL = ψ.C[0] * δˡ(0) * dag(prime(ψ.C[0])) for n in 2:dₕ temp_Rs = apply_right_transfer_matrix( @@ -304,7 +303,7 @@ function right_environment(H::InfiniteMPOMatrix, ψ::InfiniteCanonicalMPS; tol=1 return CelledVector(Rs), eᵣ[1] end -function vumps(H::InfiniteMPOMatrix, ψ::InfiniteMPS; kwargs...) +function vumps(H::InfiniteBlockMPO, ψ::InfiniteMPS; kwargs...) return vumps(H, orthogonalize(ψ, :); kwargs...) end @@ -348,7 +347,7 @@ end function tdvp_iteration_sequential( solver::Function, - H::InfiniteMPOMatrix, + H::InfiniteBlockMPO, ψ::InfiniteCanonicalMPS; (ϵᴸ!)=fill(1e-15, nsites(ψ)), (ϵᴿ!)=fill(1e-15, nsites(ψ)), @@ -410,7 +409,7 @@ end function tdvp_iteration_parallel( solver::Function, - H::InfiniteMPOMatrix, + H::InfiniteBlockMPO, ψ::InfiniteCanonicalMPS; (ϵᴸ!)=fill(1e-15, nsites(ψ)), (ϵᴿ!)=fill(1e-15, nsites(ψ)), diff --git a/test/test_iMPOConversions.jl b/test/test_iMPOConversions.jl index 2c73e83f..5f4ac39a 100644 --- a/test/test_iMPOConversions.jl +++ b/test/test_iMPOConversions.jl @@ -1,6 +1,35 @@ using ITensors using ITensorInfiniteMPS using Test + +#With the new definition of InfiniteBlockMPO, the MPO is better behaved, and hence we need to be a bit more careful +function expect_over_unit_cell(ψ::InfiniteCanonicalMPS, h::InfiniteSum{MPO}) + s = siteinds(ψ) + Ncell = nsites(h) + + energy = expect(ψ, h[1]) + for j in 2:Ncell + hf = MPO(Ncell) + for x in 1:(j - 1) + hf[x] = op("Id", s[x]) + end + for x in j:min(j + length(h[j]) - 1, Ncell) + hf[x] = h[j][x + 1 - j] + end + if j + length(h[j]) - 1 > Ncell + right_link = commonind(h[j][Ncell + 1 - j], h[j][Ncell + 2 - j]) + #dim_right = dim(right_link) + hf[end] *= onehot(dag(right_link) => 1) + elseif Ncell > j + length(h[j]) - 1 + for x in (j + length(h[j])):Ncell + hf[x] = op("Id", s[x]) + end + end + energy += expect(ψ, hf) + end + return energy +end + # # InfiniteMPO has dangling links at the end of the chain. We contract these on the outside # with l,r terminating vectors, to make a finite lattice MPO. @@ -32,6 +61,42 @@ function ITensors.expect(ψ::InfiniteCanonicalMPS, h::InfiniteMPO) return expect(ψ, terminate(h)) #defer to src/infinitecanonicalmps.jl end +function generate_edges(h::InfiniteBlockMPO) + Ncell = nsites(h) + # left termination vector + Ls = Matrix{ITensor}(undef, 1, size(h[1], 1)) + Ls[1] = ITensor(0) + for x in 2:(size(h[0], 2) - 1) + il0 = commonind(h[0][1, x], h[1][x, 1]) + Ls[x] = ITensor(0, il0) + end + Ls[end] = ITensor(1) + + Rs = Vector{ITensor}(undef, size(h[Ncell], 2)) + Rs[1] = ITensor(1) + for x in 2:(size(h[Ncell + 1], 1) - 1) + ir0 = commonind(h[Ncell + 1][x, 1], h[Ncell][1, x]) + Rs[x] = ITensor(0, ir0) + end + Rs[end] = ITensor(0) + return Ls, Rs +end + +function ITensors.expect(ψ::InfiniteCanonicalMPS, h::InfiniteBlockMPO) + Ncell = nsites(h) + L, R = generate_edges(h) + l = commoninds(ψ.AL[0], ψ.AL[1]) + L = map(a -> contract(a, δ(l, dag(prime(l)))), L) + r = commoninds(ψ.AR[Ncell + 1], ψ.AR[Ncell]) + R = map(a -> contract(a, δ(r, dag(prime(r)))), R) + L = map(a -> contract(a, ψ.C[0], dag(prime(ψ.C[0]))), L) + for j in 1:nsites(ψ) + temp = map(a -> contract(a, ψ.AR[j], dag(prime(ψ.AR[j]))), h[j]) + L = L * temp + end + return (L * R)[][] +end + #H = ΣⱼΣn (½ S⁺ⱼS⁻ⱼ₊n + ½ S⁻ⱼS⁺ⱼ₊n + SᶻⱼSᶻⱼ₊n) function ITensorInfiniteMPS.unit_cell_terms(::Model"heisenbergNNN"; NNN::Int64) opsum = OpSum() @@ -61,17 +126,6 @@ function ITensorInfiniteMPS.unit_cell_terms(::Model"hubbardNNN"; NNN::Int64) return opsum end -function ITensors.space(::SiteType"FermionK", pos::Int; p=1, q=1, conserve_momentum=true) - if !conserve_momentum - return [QN("Nf", -p) => 1, QN("Nf", q - p) => 1] - else - return [ - QN(("Nf", -p), ("NfMom", -p * pos)) => 1, - QN(("Nf", q - p), ("NfMom", (q - p) * pos)) => 1, - ] - end -end - function fermion_momentum_translator(i::Index, n::Integer; N=6) #@show n ts = tags(i) @@ -85,11 +139,7 @@ function fermion_momentum_translator(i::Index, n::Integer; N=6) return new_i end -function ITensors.op!(Op::ITensor, opname::OpName, ::SiteType"FermionK", s::Index...) - return ITensors.op!(Op, opname, SiteType("Fermion"), s...) -end - -@testset verbose = true "InfiniteMPOMatrix -> InfiniteMPO" begin +@testset verbose = true "InfiniteBlockMPO -> InfiniteMPO" begin ferro(n) = "↑" antiferro(n) = isodd(n) ? "↑" : "↓" @@ -109,11 +159,38 @@ end s = infsiteinds(site, Ncell; initstate, conserve_qns=qns) ψ = InfMPS(s, initstate) Hi = InfiniteMPO(model, s; model_kwargs...) + L = ITensor(only(commoninds(Hi[0], Hi[1]))) + R = ITensor(only(commoninds(Hi[Ncell + 1], Hi[Ncell]))) + Hi, L, R = ITensorInfiniteMPS.combineblocks_linkinds(Hi, L, R) + temp = order(L * Hi[1]) + @test temp ≈ 3 + temp = order(Hi[Ncell] * R) + @test temp ≈ 3 + Hm = InfiniteBlockMPO(model, s; model_kwargs...) + L = Array{ITensor}(undef, 1, size(Hm[1], 1)) + L[1] = ITensor() + L[end] = ITensor(1.0) + R = Vector{ITensor}(undef, size(Hm[Ncell], 2)) + R[1] = ITensor(1) + R[end] = ITensor() + for x in 2:(length(L) - 1) + L[x] = ITensor(only(commoninds(Hm[0][x, x], Hm[1][x, x]))) + R[x] = ITensor(only(commoninds(Hm[Ncell + 1][x, x], Hm[Ncell][x, x]))) + end + Hm, L, R = ITensorInfiniteMPS.combineblocks_linkinds(Hm, L, R) + temp = order.(L * Hm[1]) + ref = order.(L) .+ 2 + @test temp ≈ ref + temp = order.(Hm[Ncell] * R) + ref = order.(R) .+ 2 + @test temp ≈ ref Hs = InfiniteSum{MPO}(model, s; model_kwargs...) - Es = expect(ψ, Hs) + Es = expect_over_unit_cell(ψ, Hs) Ei = expect(ψ, Hi) + Em = expect(ψ, Hm) #@show Es Ei - @test sum(Es[1:(Ncell - NNN)]) ≈ Ei atol = 1e-14 + @test Es ≈ Ei atol = 1e-14 + @test Em ≈ Ei atol = 1e-14 end @testset "FQHE Hamitonian" begin @@ -132,10 +209,10 @@ end ψ = InfMPS(s, initstate) Hs = InfiniteSum{MPO}(model, s; model_params...) Hi = InfiniteMPO(model, s, trf; model_params...) - Es = expect(ψ, Hs) + Es = expect_over_unit_cell(ψ, Hs) Ei = expect(ψ, Hi) #@show Es Ei - @test Es[1] ≈ Ei atol = 1e-14 + @test Es ≈ Ei atol = 1e-14 end end nothing diff --git a/test/test_vumpsmpo.jl b/test/test_vumpsmpo.jl index ea2d705e..501790b2 100644 --- a/test/test_vumpsmpo.jl +++ b/test/test_vumpsmpo.jl @@ -46,7 +46,7 @@ end "sequential", "parallel" ], conserve_qns in [true, false], - nsites in [1, 2], + N in [1, 2], time_step in [-Inf] vumps_kwargs = ( @@ -58,10 +58,10 @@ end ) subspace_expansion_kwargs = (cutoff=cutoff, maxdim=maxdim) - s = infsiteinds("S=1/2", nsites; initstate, conserve_szparity=conserve_qns) + s = infsiteinds("S=1/2", N; initstate, conserve_szparity=conserve_qns) ψ = InfMPS(s, initstate) - Hmpo = InfiniteMPOMatrix(model, s; model_kwargs...) + Hmpo = InfiniteBlockMPO(model, s; model_kwargs...) # Alternate steps of running VUMPS and increasing the bond dimension ψ = tdvp(Hmpo, ψ; vumps_kwargs...) for _ in 1:outer_iters @@ -71,17 +71,15 @@ end ψ = @time tdvp(Hmpo, ψ; vumps_kwargs...) end - @test norm( - contract(ψ.AL[1:nsites]..., ψ.C[nsites]) - contract(ψ.C[0], ψ.AR[1:nsites]...) - ) ≈ 0 atol = 1e-5 - #@test contract(ψ.AL[1:nsites]..., ψ.C[nsites]) ≈ contract(ψ.C[0], ψ.AR[1:nsites]...) + @test norm(contract(ψ.AL[1:N]..., ψ.C[N]) - contract(ψ.C[0], ψ.AR[1:N]...)) ≈ 0 atol = + 1e-5 H = InfiniteSum{MPO}(model, s; model_kwargs...) energy_infinite = expect(ψ, H) - Szs_infinite = [expect(ψ, "Sz", n) for n in 1:nsites] + Szs_infinite = [expect(ψ, "Sz", n) for n in 1:N] - @test energy_finite ≈ sum(energy_infinite) / nsites rtol = 1e-4 - @test Szs_finite[nfinite:(nfinite + nsites - 1)] ≈ Szs_infinite rtol = 1e-3 + @test energy_finite ≈ sum(energy_infinite) / N rtol = 1e-4 + @test Szs_finite[nfinite:(nfinite + N - 1)] ≈ Szs_infinite rtol = 1e-3 end end @@ -135,7 +133,7 @@ end s = infsiteinds("S=1/2", nsites; conserve_szparity=conserve_qns, initstate) ψ = InfMPS(s, initstate) - Hmpo = InfiniteMPOMatrix(model, s; model_kwargs...) + Hmpo = InfiniteBlockMPO(model, s; model_kwargs...) # Alternate steps of running VUMPS and increasing the bond dimension ψ = tdvp(Hmpo, ψ; vumps_kwargs...) for _ in 1:outer_iters @@ -220,7 +218,7 @@ end ) ψ = InfMPS(s, initstate) - Hmpo = InfiniteMPOMatrix(model, s; model_kwargs...) + Hmpo = InfiniteBlockMPO(model, s; model_kwargs...) # Alternate steps of running VUMPS and increasing the bond dimension ψ = tdvp(Hmpo, ψ; vumps_kwargs...) for _ in 1:outer_iters diff --git a/test/test_vumpsmpo_fqhe.jl b/test/test_vumpsmpo_fqhe.jl index 98a156fe..e9406352 100644 --- a/test/test_vumpsmpo_fqhe.jl +++ b/test/test_vumpsmpo_fqhe.jl @@ -84,7 +84,7 @@ end ) ψ = InfMPS(s, initstate) - Hmpo = InfiniteMPOMatrix(model, s; model_kwargs...) + Hmpo = InfiniteBlockMPO(model, s; model_kwargs...) # Alternate steps of running VUMPS and increasing the bond dimension ψ = tdvp(Hmpo, ψ; vumps_kwargs...) for _ in 1:outer_iters