maartenvd / mpskit.jl Goto Github PK

Hi Maarten,

I have a question regarding how to best make use of multithreading. I read the section in the documentation which helped, but I think it does not cover my specific case.
The bulk of my simulation is time-evolution of a MPSComoving using 1-site TDVP. As far as I understand, there is nothing parallelizable at the "high-level" in MPSKit (or in general for 1-site TDVP) for this case. So the only improvement could come from low-level BLAS operations. Hence, I start Julia with only 1 julia thread, and instead set OPENBLAS_NUM_THREADS to my max number of cores (which is pretty high on our cluster = 128). This does not seem however to lead to an improvement, so I wonder what is the proper way of making use of multithreading in this case?
I suppose it boils down to how TensorKit, KrylovKit, TensorOperations handle it.
Most of them mention the Strided package.
So is the better strategy to start Julia --threads 128 and add

Strided.enable_threads()
Strided.enable_threaded_mul()

and set BLAS threads to 1?

Thanks,
Ron

leading_boundary(ψ::MPSMultiline, O::MPOMultiline, alg::GradientGrassmann) currently fails for non-trivial multiline MPOs. It seems the cost function definition for the multiline stat mech case is wrong, but it should be fairly straightforward to fix.

For example, when running

using TensorKit, MPSKit, MPSKitModels

# 2D lassical Ising
ψ = InfiniteMPS(ℂ^2, ℂ^10)
O = classical_ising(ComplexF64, Trivial; beta=1.0)

# but multiline
ψmulti = MPSMultiline([ψ, ψ])
Omulti = MPOMultiline([O, O])

alg = GradientGrassmann(; tol=1e-12)
ψ, = leading_boundary(ψmulti, Omulti, alg);

it can go wrong in any of the following ways depending on initialization:

ERROR: DomainError with -5.928428864807194e-17:
sqrt will only return a complex result if called with a complex argument. Try sqrt(Complex(x)).
Stacktrace:
 [1] throw_complex_domainerror(f::Symbol, x::Float64)
   @ Base.Math ./math.jl:33
 [2] sqrt(x::Float64)
   @ Base.Math ./math.jl:677
 [3] optimize(fg::typeof(MPSKit.GrassmannMPS.fg), x::MPSKit.GrassmannMPS.ManifoldPoint{MPSMultiline{InfiniteMPS{TrivialTensorMap{ComplexSpace, 2, 1, Matrix{ComplexF64}}, TrivialTensorMap{ComplexSpace, 1, 1, Matrix{ComplexF64}}}}, MPSKit.PerMPOInfEnv{MPOMultiline{DenseMPO{TrivialTensorMap{ComplexSpace, 2, 2, Matrix{ComplexF64}}}}, TrivialTensorMap{ComplexSpace, 2, 1, Matrix{ComplexF64}}, MPSMultiline{InfiniteMPS{TrivialTensorMap{ComplexSpace, 2, 1, Matrix{ComplexF64}}, TrivialTensorMap{ComplexSpace, 1, 1, Matrix{ComplexF64}}}}, KrylovKit.Arnoldi{KrylovKit.ModifiedGramSchmidt2, Float64}}, PeriodicArray{TensorKitManifolds.Grassmann.GrassmannTangent{TrivialTensorMap{ComplexSpace, 2, 1, Matrix{ComplexF64}}, TrivialTensorMap{ComplexSpace, 2, 1, Matrix{ComplexF64}}, TrivialTensorMap{ComplexSpace, 1, 1, Matrix{Float64}}, TrivialTensorMap{ComplexSpace, 1, 1, Matrix{ComplexF64}}}, 2}, Matrix{TrivialTensorMap{ComplexSpace, 1, 1, Matrix{ComplexF64}}}}, alg::OptimKit.ConjugateGradient{OptimKit.HagerZhang{Rational{Int64}}, Float64, OptimKit.HagerZhangLineSearch{Rational{Int64}}}; precondition::typeof(OptimKit._precondition), finalize!::typeof(OptimKit._finalize!), retract::Function, inner::typeof(MPSKit.GrassmannMPS.inner), transport!::typeof(MPSKit.GrassmannMPS.transport!), scale!::typeof(MPSKit.GrassmannMPS.scale!), add!::typeof(MPSKit.GrassmannMPS.add!), isometrictransport::Bool)
   @ OptimKit ~/.julia/packages/OptimKit/xpmbV/src/cg.jl:74

or get stuck indefinitely on

┌ Warning: mpo is not hermitian

or get stuck indefinitely on

┌ Warning: Linesearch bracket converged to a point without satisfying Wolfe conditions?
└ @ OptimKit ~/.julia/packages/OptimKit/xpmbV/src/linesearches.jl:189

Entanglementplot fails for symmetric MPSs constructed from UFCs with complex F symbols. This is because the real singular values are stored in a Vector{ComplexF64}, after which they are sorted, which fails.

using MPSKit,TensorKit,CategoryData,Plots
O = Object{UFC{3,1,2,1,1}} # Z_3
P = Vect[O](2=>1)
V = Vect[O](sector => 2 for sector in (1,2,3))
psi = InfiniteMPS(P,V)
entanglementplot(psi)

This is a follow up on #39 (comment)_ but I thought warranted a separate thread.

Changing tolgauge in TDVP, which sets tol in Lancos does not seem to have any effect.
Per the discussion in Jutho/KrylovKit.jl#52 (comment), changing eager to true leads to the expected behavior.

vumps/tdvp is parallelizeable, but @threads causes crashes every once in a while (presumeably because the whole environment-cache construction is not thread safe).

At line 105 in the finitemps.jl, I think the following code might cause some unexpected result:
(site_tensors[i],C) = leftorth!(site_tensors[i],alg=QRpos());

I think after the operation for i=1, all the other tensors are changed at the same time, which should not so. For example, consider this example:

d = 4; N = 2
hd = Int(sqrt(d))
Id = Matrix(1.0*I,hd,hd)/sqrt(2)
Id = reshape(Id, (1, d, 1))
site_tensors = [TensorMap(Id, ℝ^1 ⊗ ℝ^d, ℝ^1) for n=1:N]

i = 1
(site_tensors[i],C) = leftorth!(site_tensors[i],alg=QRpos()); normalize!(C);

In this example, site_tensors[2] is also changed! I think leftorth should replace leftorth! . Also, there may exist similar situations in other places in this package.

A possible constructor for a FiniteMPS is to specify a (vector of) physical spaces and a single maxvirtspace representing the maximum space one wants on the virtual level, according to

FiniteMPS(rand, ComplexF64, rand(3:2:20) ,pspace, maxvirtspace)

If one however wants the physical space pspace=Rep[SU₂](1//2 => 1) then this constructor will fail because it will try to find the infimum between a space containing integer reps of SU2 and another containing half-integers only. The virtual spaces of the mps should be alternating between integer and half-integer and this can currently only be done by using the constructor that requires all the virtual spaces.

when changing the data in a finitemps/finitempo/mpscomoving we should modify centerpos...

Currently, the use of virtualspace across the package can be a bit confusing, and there are contradicting choices that have been made which further add to the confusion.

The current state of affairs:

• Creating FiniteMPS objects, the constructor can take in a vector of virtual spaces which is one shorter than the length of the physical spaces, with optional kwargs for a left and right starting space. Thus, the ith element of the input vector corresponds to the virtual space to the right of site i.
• Creating InfiniteMPS objects, there is no trivial left and right starting space, thus the vector of virtual spaces is the same length as the vector fo physical spaces. In this case, the ith element of the vector corresponds to the virtual space to the left of site i.
• Querying the virtual space of an MPS object, we have left_virtualspace(state, i) and right_virtualspace(state, i), which return the left (right) virtual space of the bond tensor to the right of site i.
• entanglement_spectrum yields the singular values of the bond matrix to the right of the physical site.

Clearly, there are some inconsistencies, as well as left_virtualspace and right_virtualspace being a bit misleading.

One solution might be to:

1. Update the constructor for the InfiniteMPS to reflect that we most often refer to virtual spaces as being to the right of a physical site.
2. Introduce a new (exported) method virtualspace, which asserts that left- and right- virtualspace to the right of a site, for given state, are equal and then returns this. This should only fail if the user updated a (finite?) MPS manually with a changed bond dimension, in which case the (unexported) internal methods left_virtualspace and right_virtualspace can be used as a more expert-user-case.
3. Reflect everything clearly in the documentation and check if there are any more cases where this is relevant. In all cases, we adopt the convention that virtualspaces are defined to the right of a site index.

The current entropy function sometimes fails by returning NaN. This occurs for example after using changebonds to increase the bond dimension of a state.
Perhaps it is better to use something like this?

function entropy(state::Union{FiniteMPS,MPSComoving,InfiniteMPS},n::Int=0)
    pv = entanglement_spectrum(state,n).^2 .+ 1e-20
    -sum(log.(pv) .* pv)
end

Hello,

I encountered a strange behavior after expanding the bond dimension of an InfiniteMPS, followed by making a new InfiniteMPS using one of its tensors. I think there is something wrong with the newly created InfiniteMPS, although most operations work properly (and return same results as the original MPS), calling environments fails to converge.

I was not able to pinpoint what goes wrong exactly, but here is a minimal example:

using MPSKit,MPSKitModels,TensorKit

bondim=10

ham = nonsym_ising_ham(J=-1.0,lambda = 0.45,spin=1//2)

gs_init = InfiniteMPS([ℂ^2],[ℂ^bondim]);

(gs,_,_) = find_groundstate(gs_init,ham,VUMPS(maxiter=400));


extra_bonds = bondim
gs_expanded,_ = changebonds(gs,ham,OptimalExpand(trscheme=truncdim(extra_bonds)))

gs_2 = InfiniteMPS(gs.AL);
gs_3 = InfiniteMPS(gs_expanded.AL); #similar if using AR, AC

environments(gs_2,ham); #works ok

environments(gs_3,ham); #calcrw failed to converge 0.16636273007464153

Thanks!

Hi,
I noticed that when running memory intensive computations, Julia eats up all available RAM and does not free it up, even though it is not assigned to any variables in scope.
This sometimes leads to julia crashing. When running on a HPC cluster, my job often gets killed because it exceeds the memory given ( I request 16 GB per core).

For example, running excitation for 200 values of momenta allocates around 60GB of ram. Sometimes it crashes during this, sometimes after (during time evolution). Also had it crash at the very end of a simulation when making a plot.

I dont think there is any problem with this package, but I am wondering if you have any tips of how to deal with this.
I am new to Julia, but it seems that this issue should be quite common in these tensor networks calculations, unless I am doing something terribly wrong.
I tried playing with GC.gc() (for example, trying to free up the RAM taken by the excitation calculation) but it does not reduce it.
While I could request more RAM on the cluster, I dont see how it will help, since it seems that julia will gladly use it all!

Currently, we have this behaviour:

julia> ψ = FiniteMPS(3, ℂ^2, ℂ^1)
3-site FiniteMPS:
┌ CL[4]: TensorMap(ℂ^1 ← ℂ^1)
├── AL[3]: TensorMap((ℂ^1 ⊗ ℂ^2) ← ℂ^1)
├── AL[2]: TensorMap((ℂ^1 ⊗ ℂ^2) ← ℂ^1)
└── AL[1]: TensorMap((ℂ^1 ⊗ ℂ^2) ← ℂ^1)


julia> ψ + ψ + ψ
3-site FiniteMPS:
┌ CL[4]: TensorMap(ℂ^1 ← ℂ^1)
├── AL[3]: TensorMap((ℂ^3 ⊗ ℂ^2) ← ℂ^1)
├── AL[2]: TensorMap((ℂ^2 ⊗ ℂ^2) ← ℂ^3)
└── AL[1]: TensorMap((ℂ^1 ⊗ ℂ^2) ← ℂ^2)

While in principle, AL[3] can maximally have (ℂ^2 ⊗ ℂ^2) ← ℂ^1).

In this case, the state can in principle even be represented by a product state, but we probably do want to allow product states to be embedded in a larger bond dimension.

Hello again!

I have implemented window expansion during time evolution of an MPSComoving (along the lines of arxiv:1207.0862,arXiv: 1207.0691v2). Every timestep, I am comparing the entanglement entropy at some bond near the edge to the ground state value and if the difference is larger than some threshold I append a ground state tensor into the window.
It seems to work well but I am looking into ways to speed it up.
In particular, right now every time I extend the MPS window, I have to recalculate the environments by calling environments(mpco,th,lenvs=gs_envs,renvs=gs_envs).
I was wondering whether it is possible to reuse the environments prior to the expansion, instead of fully recalculating them?

Thanks!

finite states should include information on how they were gauged, because now time evolution algorithms assume right orthogonalized states which may be incorrect.

I'm running some GS searches for a model with 4 site unit cell and find an increasing energy as the bond dimension increases.

The context in which the 4 site unit cell arises should be irrelevant, but for completeness, it's a lattice model defined on a thin (2 site wide) strip which has unit cell of 2 in the extended direction. After snaking this becomes an effective linear model with unit cell 4.

Furthermore this behavior seems to be independent of the algorithm used to find the GS. Finally, the high D results seem to converge (from below) to some value that is not correct (I can analytically obtain the GS energy for certain parameters).

Has anyone experienced similar behavior and found a way out ?

In "function MPSKit.MPOHamiltonian(opps::SumOfLocalOperators)",
at Line 114 the variable unitcell appears to be undefined.
This creates a problem for cylinder calculations, for example :

nearest_neighbours(InfiniteCylinder(4))

ERROR: LoadError: UndefVarError: `unitcell` not defined
Stacktrace:
 [1] MPOHamiltonian(opps::SumOfLocalOperators{NTuple{8, LocalOperator{2, TensorMap{ComplexSpace, 2, 2, Trivial, Matrix{ComplexF64}, Nothing, Nothing}, LatticePoint{2, InfiniteCylinder}}}})
   @ MPSKitModels ~/work/tensorkit/MPSKitModels.jl/src/operators/mpoham.jl:114

I believe it's proper to replace Line 114 with:

if length(lattice(opps)) > 1 && _is_free_channel(data, site, hit)

When computing quasi-particle excitations a krylov algorithm is used to solve the eigenvalue problem.

According to the documentation of KrylovKit it is therefore quite risky to ask for multiple excitatations per momentum in the case that the spectrum is degenerate.

Degenerate eigenvalues

From a theoretical point of view, Krylov methods can at most find a single eigenvector associated with a targetted eigenvalue, even if the latter is degenerate. In the case of a degenerate eigenvalue, the specific eigenvector that is returned is determined by the starting vector x₀. For large problems, this turns out to be less of an issue in practice, as often a second linearly independent eigenvector is generated out of the numerical noise resulting from the orthogonalisation steps in the Lanczos or Arnoldi iteration. Nonetheless, it is important to take this into account and to try not to depend on this potentially fragile behaviour, especially for smaller problems.

As far as I know this is not mentioned in the documentation of MPSKit and on top of this it might be useful to use @warn to inform the user about this.

A more practical question : If I'd be asking for the two lowest quasiparticles which I suspect to have degenerate eigenvalues. How would this affect performance ? Could this be incredibly slow ? What exactly is ment with

generated out of numerical noise ?

and when is this expected to fail ? What happens if it does ?

It will be nice to have a push! method for SparseMPO and MPOHamiltonian so as to be able to extend the Hamiltonian when working with nonuniform Hamiltonians and a growing mps window.

Hi,

Is it possible to create and use a non uniform Hamiltonian?

Due to numerical errors, it seems like the inner product of a grassman gradient with itself can sometimes take on negative values that are of the order -eps, which then fails when OptimKit wants to take the square root.

MBA:

using MPSKit, MPSKitModels
psi = InfiniteMPS(2, 12)
H = transverse_field_ising(; hx=0.0) # exact product state solution

julia> psi, envs, delta = find_groundstate(psi, H)
[ Info: vumps @iteration 1 galerkin = 0.5339115698990319
[ Info: vumps @iteration 2 galerkin = 0.0012648372627602012
[ Info: vumps @iteration 3 galerkin = 1.1562975109758041e-11
ERROR: DomainError with -1.413494684364987e-21:
sqrt will only return a complex result if called with a complex argument. Try sqrt(Complex(x)).
Stacktrace:
 [1] throw_complex_domainerror(f::Symbol, x::Float64)
   @ Base.Math ./math.jl:33
 [2] sqrt
   @ ./math.jl:591 [inlined]
 [3] optimize(fg::typeof(MPSKit.GrassmannMPS.fg), x::MPSKit.GrassmannMPS.ManifoldPoint{InfiniteMPS{TensorKit.TensorMap{TensorKit.ComplexSpace, 2, 1, TensorKit.Trivial, Matrix{ComplexF64}, Nothing, Nothing}, TensorKit.TensorMap{TensorKit.ComplexSpace, 1, 1, TensorKit.Trivial, Matrix{ComplexF64}, Nothing, Nothing}}, MPSKit.MPOHamInfEnv{MPOHamiltonian{TensorKit.ComplexSpace, TensorKit.TensorMap{TensorKit.ComplexSpace, 2, 2, TensorKit.Trivial, Matrix{ComplexF64}, Nothing, Nothing}, ComplexF64}, TensorKit.TensorMap{TensorKit.ComplexSpace, 2, 1, TensorKit.Trivial, Matrix{ComplexF64}, Nothing, Nothing}, InfiniteMPS{TensorKit.TensorMap{TensorKit.ComplexSpace, 2, 1, TensorKit.Trivial, Matrix{ComplexF64}, Nothing, Nothing}, TensorKit.TensorMap{TensorKit.ComplexSpace, 1, 1, TensorKit.Trivial, Matrix{ComplexF64}, Nothing, Nothing}}, KrylovKit.GMRES{KrylovKit.ModifiedGramSchmidt2, Float64}}, PeriodicArray{TensorKitManifolds.Grassmann.GrassmannTangent{TensorKit.TensorMap{TensorKit.ComplexSpace, 2, 1, TensorKit.Trivial, Matrix{ComplexF64}, Nothing, Nothing}, TensorKit.TensorMap{TensorKit.ComplexSpace, 2, 1, TensorKit.Trivial, Matrix{ComplexF64}, Nothing, Nothing}, TensorKit.TensorMap{TensorKit.ComplexSpace, 1, 1, TensorKit.Trivial, Matrix{Float64}, Nothing, Nothing}, TensorKit.TensorMap{TensorKit.ComplexSpace, 1, 1, TensorKit.Trivial, Matrix{ComplexF64}, Nothing, Nothing}}, 1}, Vector{TensorKit.TensorMap{TensorKit.ComplexSpace, 1, 1, TensorKit.Trivial, Matrix{ComplexF64}, Nothing, Nothing}}}, alg::OptimKit.ConjugateGradient{OptimKit.HagerZhang{Rational{Int64}}, Float64, OptimKit.HagerZhangLineSearch{Rational{Int64}}}; precondition::typeof(OptimKit._precondition), finalize!::typeof(OptimKit._finalize!), retract::Function, inner::typeof(MPSKit.GrassmannMPS.inner), transport!::typeof(MPSKit.GrassmannMPS.transport!), scale!::typeof(MPSKit.GrassmannMPS.scale!), add!::typeof(MPSKit.GrassmannMPS.add!), isometrictransport::Bool)
   @ OptimKit ~/.julia/packages/OptimKit/xpmbV/src/cg.jl:27
 [4] find_groundstate(Ψ::InfiniteMPS{TensorKit.TensorMap{TensorKit.ComplexSpace, 2, 1, TensorKit.Trivial, Matrix{ComplexF64}, Nothing, Nothing}, TensorKit.TensorMap{TensorKit.ComplexSpace, 1, 1, TensorKit.Trivial, Matrix{ComplexF64}, Nothing, Nothing}}, H::MPOHamiltonian{TensorKit.ComplexSpace, TensorKit.TensorMap{TensorKit.ComplexSpace, 2, 2, TensorKit.Trivial, Matrix{ComplexF64}, Nothing, Nothing}, ComplexF64}, alg::GradientGrassmann, envs::MPSKit.MPOHamInfEnv{MPOHamiltonian{TensorKit.ComplexSpace, TensorKit.TensorMap{TensorKit.ComplexSpace, 2, 2, TensorKit.Trivial, Matrix{ComplexF64}, Nothing, Nothing}, ComplexF64}, TensorKit.TensorMap{TensorKit.ComplexSpace, 2, 1, TensorKit.Trivial, Matrix{ComplexF64}, Nothing, Nothing}, InfiniteMPS{TensorKit.TensorMap{TensorKit.ComplexSpace, 2, 1, TensorKit.Trivial, Matrix{ComplexF64}, Nothing, Nothing}, TensorKit.TensorMap{TensorKit.ComplexSpace, 1, 1, TensorKit.Trivial, Matrix{ComplexF64}, Nothing, Nothing}}, KrylovKit.GMRES{KrylovKit.ModifiedGramSchmidt2, Float64}})
   @ MPSKit ~/Projects/JuliaProjects/MPSKit.jl/src/algorithms/groundstate/gradient_grassmann.jl:53
 [5] find_groundstate
   @ ~/Projects/JuliaProjects/MPSKit.jl/src/algorithms/unionalg.jl:25 [inlined]
 [6] #find_groundstate#381
   @ ~/Projects/JuliaProjects/MPSKit.jl/src/algorithms/groundstate/find_groundstate.jl:41 [inlined]
 [7] find_groundstate (repeats 2 times)
   @ ~/Projects/JuliaProjects/MPSKit.jl/src/algorithms/groundstate/find_groundstate.jl:19 [inlined]
 [8] top-level scope
   @ REPL[12]:1

It might be useful to add a dedicated approximate(Ψ₀::AbstractMPS, Ψ::AbstractMPS, algorithm, [environments]; kwargs...) method for a globally optimal truncation of a given MPS. It is mentioned in the changebonds documentation that combining changebonds(Ψ::AbstractMPS, alg::SvdCt) with approximate can give a globally optimal truncation, but as far as I can tell all current implementations of approximate take an MPO-MPS pair as the object to approximate, so I don't see how this could be done directly.

Currently the default algorithm for singular value decomposition is a divide and conquer algorithm.

From an optimization perspective this is great. Indeed, for large tensors such algorithms outperform the standard method Householder(?) method but they are known to be less stable and prone to crashing.

One specific stack trace for such a crash in MPSKit might be :

but I suspect there are a multiple a ways to reproduce this crash.

In my particular case the issue was patched by changing projectisometric! from TensorKitManifolds.jl into :

function projectisometric!(W::AbstractTensorMap; alg=Polar())
    if alg isa TensorKit.Polar || alg isa TensorKit.SDD
        foreach(blocks(W)) do (c, b)
            #sometimes the divide and conquer algorithm used in polarsdd fails and gives LAPACK error 1
            #in that case, we try the standard SVD algorithm which is slightly more expensive but more robust
            #in case that fails I'll be a sad camper... (spoiler, it doesn't)
            b_backup = copy(b)
            try 
                return _polarsdd!(b)
            catch
                b = b_backup
                return _polarsvd!(b)
            end
        end
    elseif alg isa TensorKit.SVD
        foreach(blocks(W)) do (c, b)
            return _polarsvd!(b)
        end
    elseif alg isa PolarNewton
        foreach(blocks(W)) do (c, b)
            return _polarnewton!(b)
        end
    else
        throw(ArgumentError("unkown algorithm for projectisometric!: alg = $alg"))
    end
    return W
end

so that the slower, more stable method is used whenever the initial attempt fails.

Clearly this solution is quite (one could even say very) messy and slow (since it has to make a backup of the to-be-svd'd tensor) so that much improvement is still possible. In particular MPSKit should offer the option to choose different SVD decompositions on default. And perhaps there is a smarter way to implement the try cath loop without the necessity of a backup ?

When running finite excitation it is possible to get the following error

MethodError: no method matching expectation_value(::FiniteMPS{TensorMap{GradedSpace{PlanarTrivial, Tuple{Int64}}, 2, 1, PlanarTrivial, TensorKit.SortedVectorDict{PlanarTrivial, Matrix{ComplexF64}}, FusionTree{PlanarTrivial, 2, 0, 1, Nothing}, FusionTree{PlanarTrivial, 1, 0, 0, Nothing}}, TensorMap{GradedSpace{PlanarTrivial, Tuple{Int64}}, 1, 1, PlanarTrivial, TensorKit.SortedVectorDict{PlanarTrivial, Matrix{ComplexF64}}, FusionTree{PlanarTrivial, 1, 0, 0, Nothing}, FusionTree{PlanarTrivial, 1, 0, 0, Nothing}}}, ::MPSKit.LinearCombination{Tuple{MPOHamiltonian{GradedSpace{PlanarTrivial, Tuple{Int64}}, TensorMap{GradedSpace{PlanarTrivial, Tuple{Int64}}, 2, 2, PlanarTrivial, TensorKit.SortedVectorDict{PlanarTrivial, Matrix{ComplexF64}}, FusionTree{PlanarTrivial, 2, 0, 1, Nothing}, FusionTree{PlanarTrivial, 2, 0, 1, Nothing}}, ComplexF64}, MPSKit.ProjectionOperator{FiniteMPS{TensorMap{GradedSpace{PlanarTrivial, Tuple{Int64}}, 2, 1, PlanarTrivial, TensorKit.SortedVectorDict{PlanarTrivial, Matrix{ComplexF64}}, FusionTree{PlanarTrivial, 2, 0, 1, Nothing}, FusionTree{PlanarTrivial, 1, 0, 0, Nothing}}, TensorMap{GradedSpace{PlanarTrivial, Tuple{Int64}}, 1, 1, PlanarTrivial, TensorKit.SortedVectorDict{PlanarTrivial, Matrix{ComplexF64}}, FusionTree{PlanarTrivial, 1, 0, 0, Nothing}, FusionTree{PlanarTrivial, 1, 0, 0, Nothing}}}}}, Tuple{Float64, Float64}}, ::MPSKit.LazyLincoCache{MPSKit.LinearCombination{Tuple{MPOHamiltonian{GradedSpace{PlanarTrivial, Tuple{Int64}}, TensorMap{GradedSpace{PlanarTrivial, Tuple{Int64}}, 2, 2, PlanarTrivial, TensorKit.SortedVectorDict{PlanarTrivial, Matrix{ComplexF64}}, FusionTree{PlanarTrivial, 2, 0, 1, Nothing}, FusionTree{PlanarTrivial, 2, 0, 1, Nothing}}, ComplexF64}, MPSKit.ProjectionOperator{FiniteMPS{TensorMap{GradedSpace{PlanarTrivial, Tuple{Int64}}, 2, 1, PlanarTrivial, TensorKit.SortedVectorDict{PlanarTrivial, Matrix{ComplexF64}}, FusionTree{PlanarTrivial, 2, 0, 1, Nothing}, FusionTree{PlanarTrivial, 1, 0, 0, Nothing}}, TensorMap{GradedSpace{PlanarTrivial, Tuple{Int64}}, 1, 1, PlanarTrivial, TensorKit.SortedVectorDict{PlanarTrivial, Matrix{ComplexF64}}, FusionTree{PlanarTrivial, 1, 0, 0, Nothing}, FusionTree{PlanarTrivial, 1, 0, 0, Nothing}}}}}, Tuple{Float64, Float64}}, Tuple{MPSKit.FinEnv{Nothing, MPOHamiltonian{GradedSpace{PlanarTrivial, Tuple{Int64}}, TensorMap{GradedSpace{PlanarTrivial, Tuple{Int64}}, 2, 2, PlanarTrivial, TensorKit.SortedVectorDict{PlanarTrivial, Matrix{ComplexF64}}, FusionTree{PlanarTrivial, 2, 0, 1, Nothing}, FusionTree{PlanarTrivial, 2, 0, 1, Nothing}}, ComplexF64}, TensorMap{GradedSpace{PlanarTrivial, Tuple{Int64}}, 2, 1, PlanarTrivial, TensorKit.SortedVectorDict{PlanarTrivial, Matrix{ComplexF64}}, FusionTree{PlanarTrivial, 2, 0, 1, Nothing}, FusionTree{PlanarTrivial, 1, 0, 0, Nothing}}, Vector{TensorMap{GradedSpace{PlanarTrivial, Tuple{Int64}}, 2, 1, PlanarTrivial, TensorKit.SortedVectorDict{PlanarTrivial, Matrix{ComplexF64}}, FusionTree{PlanarTrivial, 2, 0, 1, Nothing}, FusionTree{PlanarTrivial, 1, 0, 0, Nothing}}}}, MPSKit.FinEnv{FiniteMPS{TensorMap{GradedSpace{PlanarTrivial, Tuple{Int64}}, 2, 1, PlanarTrivial, TensorKit.SortedVectorDict{PlanarTrivial, Matrix{ComplexF64}}, FusionTree{PlanarTrivial, 2, 0, 1, Nothing}, FusionTree{PlanarTrivial, 1, 0, 0, Nothing}}, TensorMap{GradedSpace{PlanarTrivial, Tuple{Int64}}, 1, 1, PlanarTrivial, TensorKit.SortedVectorDict{PlanarTrivial, Matrix{ComplexF64}}, FusionTree{PlanarTrivial, 1, 0, 0, Nothing}, FusionTree{PlanarTrivial, 1, 0, 0, Nothing}}}, Vector{Nothing}, TensorMap{GradedSpace{PlanarTrivial, Tuple{Int64}}, 2, 1, PlanarTrivial, TensorKit.SortedVectorDict{PlanarTrivial, Matrix{ComplexF64}}, FusionTree{PlanarTrivial, 2, 0, 1, Nothing}, FusionTree{PlanarTrivial, 1, 0, 0, Nothing}}, TensorMap{GradedSpace{PlanarTrivial, Tuple{Int64}}, 1, 3, PlanarTrivial, TensorKit.SortedVectorDict{PlanarTrivial, Matrix{ComplexF64}}, FusionTree{PlanarTrivial, 1, 0, 0, Nothing}, FusionTree{PlanarTrivial, 3, 1, 2, Nothing}}}}})

The error details that there does not exist a method for expectation_value(::FiniteMPS,::MPSKit.LinearCombination,::MPSKit.LazyLincoCache). More precisely the error stems from this part inside dmrgexciation.jl

super_op = LinearCombination(
        tuple(H, ProjectionOperator.(states)...),
        tuple(1.0, broadcast(x -> alg.weight, states)...),
    )
    envs = environments(init, super_op)
    ne, _ = find_groundstate(init, super_op, alg.gsalg, envs)

When the maximum amount of iterations is reached find_groundstate warns the users with a message that contains expectation_value of that state with H. In this case H= super_op is a LinearCombination so it tries to calculate this but fails because this method is not defined. A possible solution is to define

function expectation_value(Ψ, H::LinearCombination, envs::LazyLincoCache=environments(Ψ,H))
    return return sum(((c, op, env),) -> c * expectation_value(Ψ, op, env), zip(H.coeffs,H.opps, envs.envs))
end

However this causes again an error when it tries to do expectation_value(Ψ, ::MPSKit.ProjectionOperator, ::MPSKit.FinEnv) for which there also doesn't exist any method. I don't know how expectation_value should be defined for a ProjectionOperator and at this point it seems easier to just suppress the warning that causes the error in the first place.

Having a look at the output of the tests, it seems like there is some performance issue going on with the combination of LazySum and GradientGrassmann. My best guess is that the gradient is actually not computed entirely correctly, but the algorithm still converges because Hlazy = [0.5*H - H + 5.553H] is actually not that good of a testcase.

Test Summary:        | Pass  Total     Time
find_groundstate     |   40     40  1m16.4s
  Infinite 1         |    2      2     0.2s
  Infinite 2         |    2      2     0.8s
  Infinite 3         |    2      2     3.2s
  Infinite 4         |    2      2     7.1s
  Infinite 5         |    2      2     0.2s
  LazySum Infinite 1 |    3      3     2.2s
  LazySum Infinite 2 |    3      3     3.1s
  LazySum Infinite 3 |    3      3     2.9s
  LazySum Infinite 4 |    3      3     5.3s
  LazySum Infinite 5 |    3      3     0.9s
  Finite 1           |    2      2     1.3s
  Finite 2           |    2      2     2.6s
  Finite 3           |    2      2     7.2s
  LazySum Finite 1   |    3      3     2.9s
  LazySum Finite 2   |    3      3     3.0s
  LazySum Finite 3   |    3      3    33.4s

Hi,

Is it possible to increase the size of the "B" tensor in the quasiparticle ansatz to more than a single site ?

My understanding is that in some cases a proper description of a quasiparticle may require more than one site, and the right way of checking for convergence (https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.111.080401) is by increasing the size of the "B" block rather than its bond dimension.

given a state with Float entries and a Complex hamiltonian, timestep! errors because it tries to update in-place. Probably could do with a better check and error message

Because correlator is written using @tensor, it produces incorrect results for non-bosonic braiding rules

I'm not sure how to convert the MPO representation as Vector of the Hamiltonian for ex. the transverse Ising chain into MPSKit's MPOHamiltonian. I think it can be casted into a periodic array but I'm not sure about the details. Are there perhaps examples? Many thanks!

Hi,
Thanks for this great package!

I noticed that the quasiparticle (QP, InfiniteQP, FiniteQP) types are not exported. Can they be exported?
It would help in writing functions that accept such objects (obtained via the excitations function).

On that subject, may I ask what is the purpose of the "utility leg" of the B tensor in the quasiparticle state?
I find myself having to manually remove it (and since I could not find an easy way to accomplish it with TensorKit, I just convert it to an array and build a new TensorMap without the useless index).

Hi,

I think there might be a problem with the gauge conversion of the quasiparticle type, when the left and right ground states are different.
The resulting B tensor does not satisfy the right gauge-fixing condition.

Here is a minimal example:

using MPSKit,MPSKitModels
using TensorKit


########### helper functions to find symmetry broken ground states
function get_bond_dimensions(state::Union{FiniteMPS,InfiniteMPS,MPSComoving}) 
    N=length(state)+1
    bonds = Array{Int}(undef,N)
    for i in 1:N
        bonds[i] = dim(left_virtualspace(state,i-1))
    end
    return bonds
end

get_max_bond_dim(state::Union{FiniteMPS,InfiniteMPS,MPSComoving}) = maximum(get_bond_dimensions(state))

function VUMPS_with_init_state(init_gs_guess::InfiniteMPS,ham::MPOHamiltonian,max_bond_dim::Int;maxiter::Int = 400)
    function my_finalize(iter,state,ham,envs)
        if get_max_bond_dim(state) < max_bond_dim::Int
            state,envs = changebonds(state,ham,VUMPSSvdCut(trscheme = truncdim(max_bond_dim)),envs)
        end
        return state,envs;
    end
    return find_groundstate(init_gs_guess,ham,VUMPS(maxiter=maxiter,finalize=my_finalize));
end 
###########

ham = nonsym_ising_ham(J=-1.0,lambda = 0.5,spin=1//2)

sx,sy,sz  =nonsym_spintensors(1//2)

up = zeros(ComplexF64,1,2,1);
up[:,1,:] .= 1.0
up_data = TensorMap(up,ℂ^1 ⊗ ℂ^2,ℂ^1);
dn = zeros(ComplexF64,1,2,1);
dn[:,2,:] .= 1.0
dn_data = TensorMap(dn,ℂ^1 ⊗ ℂ^2,ℂ^1);
up_state = InfiniteMPS([up_data]);
dn_state = InfiniteMPS([dn_data]);

# find the two symmetry broken ground states
gs⁺,_ = VUMPS_with_init_state(up_state,ham,10);
gs⁻,_ = VUMPS_with_init_state(dn_state,ham,10);

# check
@info expectation_value(gs⁺,sz)
@info expectation_value(gs⁻,sz)

# find kink excitation
k = 0.2;
#(triv_energies,triv_Bs) = excitations(ham,QuasiparticleAnsatz(),k,gs⁺,num=1);
(top_energies_pm,top_Bs_pm) = excitations(ham,QuasiparticleAnsatz(),k,gs⁺,environments(gs⁺,ham),gs⁻,environments(gs⁻,ham),num=1);


qp_state = top_Bs_pm[1];
qp_stateR = convert(RightGaugedQP,top_Bs_pm[1]);


B = qp_state[1];
BR = qp_stateR[1];

AL = qp_state.left_gs.AL[1]
AR = qp_stateR.right_gs.AR[1]

@tensor LGzero[k l; m] := B[i j; k l]*conj(AL[i j;m]);

@info norm(LGzero) # 2.8262689435084304e-16


@tensor RGzero[k l; t] := BR[l j; k m]*conj(AR[t j;m]);

@info norm(RGzero) #0.24027386103287698

On

Hi Maarten,

would you be open to some suggestions for renaming things.

• all names with Mps and Mpo to MPS and MPO
• MPSCenterGauged to just InfiniteMPS, the details of how this is stored/represented are not really relevant to most users I assume
• Is it possible to fuse FiniteMPS and MPSComoving? The ground state and time evolution algorithms anyway accept Union{FiniteMPS,MPSComoving}

Is there an advantage to having FiniteMPS <: AbstractArray?

I can prepare a PR somewhere next week.

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

Because the setup.jl file is included for each of the test groups, but none of these actually use modules to separate their namespaces, julia throws a bunch of complaints about re-including a file. This could easily be fixed by making either/both setup.jl and the test groups their own module.

When performing groundstate searches requiring many iterations it is desirable to have some notion of convergence speed to estimate runtime on clusters and better help planning.

I used to have an implementation of this using UnicodePlots.jl which can run in any terminal.

I was thinking about re-implementing this and was wondering how such functionality would fit most elegantly within the current framework. I suspect the finalize function can be used for this ? Alternatively it might be implemented as a verbosity setting.

When starting with an InfiniteMPS from tensors that are not full rank, the gauging algorithm does not necessarily lead to square CR tensors. If this happens at the edges of the unitcell, IDMRG2 fails because it takes an inverse of C which is not square.

julia> using MPSKit, TensorKit
julia> As = [MPSTensor(ℂ^2, ℂ^1, ℂ^3), MPSTensor(ℂ^2, ℂ^3, ℂ^4), MPSTensor(ℂ^2, ℂ^4, ℂ^1)];

julia> psi = InfiniteMPS(As)
3-site InfiniteMPS:
│   ⋮
│ CR[3]: TensorMap(ℂ^1 ← ℂ^1)
├── AL[3]: TensorMap((ℂ^4 ⊗ ℂ^2) ← ℂ^1)
├── AL[2]: TensorMap((ℂ^2 ⊗ ℂ^2) ← ℂ^4)
├── AL[1]: TensorMap((ℂ^1 ⊗ ℂ^2) ← ℂ^2)
│   ⋮


julia> space.(psi.CR)
3-element PeriodicArray{TensorMapSpace{ComplexSpace, 1, 1}, 1}:
 ℂ^2 ← ℂ^2
 ℂ^4 ← ℂ^2
 ℂ^1 ← ℂ^1

In general, there are many subtle problems for MPS tensors that aren't full rank, so I think we should both add an explicit check in the InfiniteMPS constructor that prohibits this, as well as provide a mechanism for creating full rank MPS tensors from non-full rank ones.

partial stack trace for completeness:

LoadError: SpaceMismatch("codomain ProductSpace(Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[SU₂])]((0, 0, 0)=>16, (0, 4, 0)=>35, (0, -4, 0)=>1, (0, 8, 0)=>16, (0, 12, 0)=>1, (1, 2, 1/2)=>44, (1, -2, 1/2)=>8, (1, 6, 1/2)=>44, (1, 10, 1/2)=>8, (0, 0, 1)=>22, (0, 4, 1)=>51, (0, -4, 1)=>1, (0, 8, 1)=>22, (0, 12, 1)=>1, (1, 2, 3/2)=>28, (1, -2, 3/2)=>4, (1, 6, 3/2)=>28, (1, 10, 3/2)=>4, (0, 0, 2)=>6, (0, 4, 2)=>17, (0, 8, 2)=>6, (1, 2, 5/2)=>4, (1, 6, 5/2)=>4, (0, 4, 3)=>1)) and domain ProductSpace(Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[SU₂])]((0, 0, 0)=>3, (0, 4, 0)=>1, (0, -4, 0)=>1, (1, 2, 1/2)=>2, (1, -2, 1/2)=>2, (0, 0, 1)=>1)) are not isomorphic: no inverse")
Stacktrace:
 [1] inv(t::TensorMap{GradedSpace{ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, TensorKit.SortedVectorDict{ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, Int64}}, 1, 1, ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, TensorKit.SortedVectorDict{ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, Matrix{ComplexF64}}, FusionTree{ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, 1, 0, 0, Nothing}, FusionTree{ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, 1, 0, 0, Nothing}})
   @ TensorKit ~/.julia/packages/TensorKit/gpuAY/src/tensors/linalg.jl:268
 [2] find_groundstate(ost::InfiniteMPS{TensorMap{GradedSpace{ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, TensorKit.SortedVectorDict{ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, Int64}}, 2, 1, ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, TensorKit.SortedVectorDict{ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, Matrix{ComplexF64}}, FusionTree{ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, 2, 0, 1, Nothing}, FusionTree{ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, 1, 0, 0, Nothing}}, TensorMap{GradedSpace{ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, TensorKit.SortedVectorDict{ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, Int64}}, 1, 1, ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, TensorKit.SortedVectorDict{ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, Matrix{ComplexF64}}, FusionTree{ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, 1, 0, 0, Nothing}, FusionTree{ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, 1, 0, 0, Nothing}}}, H::MPOHamiltonian{GradedSpace{ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, TensorKit.SortedVectorDict{ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, Int64}}, TensorMap{GradedSpace{ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, TensorKit.SortedVectorDict{ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, Int64}}, 2, 2, ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, TensorKit.SortedVectorDict{ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, Matrix{ComplexF64}}, FusionTree{ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, 2, 0, 1, Nothing}, FusionTree{ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, 2, 0, 1, Nothing}}, ComplexF64}, alg::IDMRG2{KrylovKit.Arnoldi{KrylovKit.ModifiedGramSchmidt2, Float64}}, oenvs::MPSKit.MPOHamInfEnv{MPOHamiltonian{GradedSpace{ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, TensorKit.SortedVectorDict{ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, Int64}}, TensorMap{GradedSpace{ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, TensorKit.SortedVectorDict{ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, Int64}}, 2, 2, ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, TensorKit.SortedVectorDict{ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, Matrix{ComplexF64}}, FusionTree{ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, 2, 0, 1, Nothing}, FusionTree{ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, 2, 0, 1, Nothing}}, ComplexF64}, TensorMap{GradedSpace{ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, TensorKit.SortedVectorDict{ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, Int64}}, 2, 1, ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, TensorKit.SortedVectorDict{ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, Matrix{ComplexF64}}, FusionTree{ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, 2, 0, 1, Nothing}, FusionTree{ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, 1, 0, 0, Nothing}}, InfiniteMPS{TensorMap{GradedSpace{ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, TensorKit.SortedVectorDict{ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, Int64}}, 2, 1, ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, TensorKit.SortedVectorDict{ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, Matrix{ComplexF64}}, FusionTree{ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, 2, 0, 1, Nothing}, FusionTree{ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, 1, 0, 0, Nothing}}, TensorMap{GradedSpace{ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, TensorKit.SortedVectorDict{ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, Int64}}, 1, 1, ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, TensorKit.SortedVectorDict{ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, Matrix{ComplexF64}}, FusionTree{ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, 1, 0, 0, Nothing}, FusionTree{ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, 1, 0, 0, Nothing}}}, KrylovKit.GMRES{KrylovKit.ModifiedGramSchmidt2, Float64}})
   @ MPSKit ~/.julia/packages/MPSKit/JT8u6/src/algorithms/groundstate/idmrg.jl:116
 [3] find_groundstate(ost::InfiniteMPS{TensorMap{GradedSpace{ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, TensorKit.SortedVectorDict{ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, Int64}}, 2, 1, ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, TensorKit.SortedVectorDict{ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, Matrix{ComplexF64}}, FusionTree{ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, 2, 0, 1, Nothing}, FusionTree{ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, 1, 0, 0, Nothing}}, TensorMap{GradedSpace{ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, TensorKit.SortedVectorDict{ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, Int64}}, 1, 1, ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, TensorKit.SortedVectorDict{ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, Matrix{ComplexF64}}, FusionTree{ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, 1, 0, 0, Nothing}, FusionTree{ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, 1, 0, 0, Nothing}}}, H::MPOHamiltonian{GradedSpace{ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, TensorKit.SortedVectorDict{ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, Int64}}, TensorMap{GradedSpace{ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, TensorKit.SortedVectorDict{ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, Int64}}, 2, 2, ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, TensorKit.SortedVectorDict{ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, Matrix{ComplexF64}}, FusionTree{ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, 2, 0, 1, Nothing}, FusionTree{ProductSector{Tuple{FermionParity, U1Irrep, SU2Irrep}}, 2, 0, 1, Nothing}}, ComplexF64}, alg::IDMRG2{KrylovKit.Arnoldi{KrylovKit.ModifiedGramSchmidt2, Float64}})
   @ MPSKit ~/.julia/packages/MPSKit/JT8u6/src/algorithms/groundstate/idmrg.jl:85

When calling environments on a WindowMPS without specifying the lenvs and renvskwargs it will calculate these with the supplied Hamiltonian. For some cases this will result in an AssertionError: len == length(ham) down the line because at some point an explicit check for MPOHamInfEnv that the length of H and the InfiniteMPS need to be equal, was introduced. For example these cases will result in an error.

Ψ = InfiniteMPS([ℂ^2],[ℂ^2]);
window = WindowMPS(Ψ,20);
nn = TensorMap(rand, ComplexF64, ℂ^2 * ℂ^2, ℂ^2 * ℂ^2);
nn += nn'
H = repeat(MPOHamiltonian(nn),20);
environments(window,H)

Ψl = InfiniteMPS([ℂ^2],[ℂ^2]);
Ψf = FiniteMPS(repeat([ℂ^2],20),repeat([ℂ^2],21))
Ψr = repeat(Ψl,2);
nn = TensorMap(rand, ComplexF64, ℂ^2 * ℂ^2, ℂ^2 * ℂ^2);
nn += nn'
H = repeat(MPOHamiltonian(nn),2);
window2 = WindowMPS(Ψl,Ψf,Ψr);
environments(window,H)

A potential solution would be to make the specification of lenv and renv mandatory or demand separate Hamiltonians for the different parts of the window.

Hi, I notice the MPO in MPSKit currently only supports single-site unit-cell models. Is it possible to construct MPOs with a 2-site unit cell? For example, there are different odd and even bond interactions in a Hamiltonian or finite width ladder systems. Thanks.

approximate(ψ::FiniteMPS, (H::DensMPO, ψ₀::FiniteMPS), alg::Union{DMRG,DMRG2}) is currently broken. It seems there's a method missing to automatically construct the leftstart and rightstart for the environment initialization in the case of a DenseMPO.

using TensorKit, MPSKit

N = 4
ψ₀ = FiniteMPS(randn, ComplexF64, fill(ℂ^2, N), [ℂ^1, fill(ℂ^5, N-1)..., ℂ^1])

O = TensorMap(randn, ComplexF64, ℂ^2 ⊗ ℂ^2 ← ℂ^2 ⊗ ℂ^2)
OL = TensorMap(randn, ComplexF64, ℂ^1 ⊗ ℂ^2 ← ℂ^2 ⊗ ℂ^2)
OR = TensorMap(randn, ComplexF64, ℂ^2 ⊗ ℂ^2 ← ℂ^2 ⊗ ℂ^1)
H = DenseMPO([OL, fill(O, N-2)..., OR])

ψ = FiniteMPS(randn, ComplexF64, fill(ℂ^2, N), [ℂ^1, fill(ℂ^5, N-1)..., ℂ^1])
ψ, = approximate(ψ, (H, ψ₀), DMRG())

throws

ERROR: MethodError: no method matching environments(::FiniteMPS{TrivialTensorMap{ComplexSpace, 2, 1, Matrix{ComplexF64}}, TrivialTensorMap{ComplexSpace, 1, 1, Matrix{ComplexF64}}}, ::DenseMPO{TrivialTensorMap{ComplexSpace, 2, 2, Matrix{ComplexF64}}}, ::FiniteMPS{TrivialTensorMap{ComplexSpace, 2, 1, Matrix{ComplexF64}}, TrivialTensorMap{ComplexSpace, 1, 1, Matrix{ComplexF64}}})

Closest candidates are:
  environments(::Any, ::Any, ::Any, ::Any)
   @ MPSKit ~/.julia/packages/MPSKit/t8tLc/src/environments/FinEnv.jl:24
  environments(::Any, ::Any, ::Any, ::Any, ::Any)
   @ MPSKit ~/.julia/packages/MPSKit/t8tLc/src/environments/FinEnv.jl:27
  environments(::FiniteMPS{S}, ::Union{MPOHamiltonian, SparseMPO}, ::Any) where S
   @ MPSKit ~/.julia/packages/MPSKit/t8tLc/src/environments/FinEnv.jl:47
  ...

I am wondering whether it is possible to write all the AL,AR,CR,AC tensors and even the environment tensors of a MPS (sometimes the unit cell can be large) into a single file? As MPSKit relies on TensorKit, I think one has to store each TensorMap into an independent file for now, which is not convenient sometimes. So is there a simple way to realize it?