Change sector convention for excitations

[lkdvos]

Change sector convention to be more intuitive, in particular allowing for:

eigvals(convert(TensorMap, H))[sector] == exact_diagonalization(H; sector)

and equivalent results for the excitations.

to do

• ED
• excitations
• add_physical_charge
• rebase on Improve quasiparticle type stability #272

Sidenote

While doing this I discovered that the excitation code was type unstable for some reasons, which I also resolved, #272 should be merged first

[codecov]

Codecov Report

Attention: Patch coverage is 98.00000% with 1 line in your changes missing coverage. Please review.

Files with missing lines	Patch %	Lines
src/states/finitemps.jl	94.11%	1 Missing ⚠️

Files with missing lines	Coverage Δ
src/algorithms/ED.jl	100.00% <100.00%> (+100.00%)	⬆️
src/operators/abstractmpo.jl	67.91% <100.00%> (+10.69%)	⬆️
src/states/quasiparticle_state.jl	87.14% <100.00%> (+0.10%)	⬆️
src/states/finitemps.jl	80.27% <94.11%> (+2.30%)	⬆️

... and 3 files with indirect coverage changes

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[leburgel]

Looks perfect to me, but it's a bit annoying that really none of the changed parts were covered by the tests to begin with so we're essentially eyeballing to check that things work. The PR itself achieves the goal though, so I would be okay with merging despite the lack of tests and add these in a follow-up.

As a final comment, it might be good to list some necessary downstream changes to account for the changes here while we're at it, just to keep track of things. Is there anything else besides MPSKitModels.jl and PEPSKit.jl?

[lkdvos]

I don't think there's anything else besides MPSKitModels and PEPSKit. Did you have anything particular in mind for listing these? It's a bit strange to open issues in these packages for things that will be broken once the breaking release of MPSKit comes out, which won't even be "breaking" since these packages shouldn't auto-update because of compat. On the other hand, probably this doesn't matter that much so we might as well do that.

I've added some more tests, which should now explicitly check at least some of these conventions (and also the exact_diagonalization function).

[leburgel]

Thanks for the effort of adding tests right away!

For listing the required downstream changes, I just meant adding some note in the PR description here, so that there's a record of what should be done and we can refer to this PR when updating the compatibility of downstream packages. My question was a bit unclear, but I didn't mean opening actual issues in the downstream repos. I also don't think that's necessary.

[ZongYongyue]

Hi Lukas,

I really like this part of the improvement, but I encountered an issue when trying to use the new method to construct a charged MPS. When I place one electron on the right side, everything works correctly—both the ground state energy and particle number match expectations. However, when I place two electrons, the results remain the same as with one electron. Could you please help me identify where the issue might be?

I tested this on a Hubbard model with size = 4. The first half of the implementation uses the new way, while the second half uses my old way for benchmarking:

using TensorKit
using MPSKit
using MPSKitModels:FiniteChain
using DynamicalCorrelators: hubbard, randFiniteMPS

N = 4
elt = Float64

I = FermionParity ⊠ SU2Irrep ⊠ U1Irrep 
pspace = Vect[I]((0, 0, -1) => 1, (1, 1//2, 0) => 1, (0, 0, 1) => 1)
right=Vect[I]((1, 1//2, 1)=> 2)
H = hubbard(Float64, SU2Irrep, U1Irrep, FiniteChain(N); t=1, U=8, μ=0)
st = FiniteMPS(MPSKit.physicalspace(H), MPSKit.max_virtualspaces(MPSKit.physicalspace(H); right)[2:(end - 1)];right)
gs, envs, delta = find_groundstate(st, H, DMRG2(trscheme= truncerr(1e-6)));
E0 = expectation_value(gs, H)
n = TensorMap(zeros, elt, pspace ← pspace)
block(n, I((0,0,1))) .= 2
block(n, I((1,1//2,0))) .= 1
id = isomorphism(n)
num = sum([correlator(gs, MPSKit.add_util_leg(n), MPSKit.add_util_leg(id), 1, 2),
        correlator(gs, MPSKit.add_util_leg(n), MPSKit.add_util_leg(id), 2, 3),
        correlator(gs, MPSKit.add_util_leg(n), MPSKit.add_util_leg(id), 3, 4),
        correlator(gs, MPSKit.add_util_leg(id), MPSKit.add_util_leg(n), 1, 4)])
println("E0: ", real(E0),"  Particle number: ", real(num))

filling = (3,2)
P, Q = filling
st = randFiniteMPS(elt, SU2Irrep, U1Irrep, N; filling=filling)
H = hubbard(elt, SU2Irrep, U1Irrep, FiniteChain(4); filling=filling, t=1, U=8, μ=0)
gs, envs, delta = find_groundstate(st, H, DMRG2(trscheme= truncerr(1e-6)));
E0 = expectation_value(gs, H)
dim.(MPSKit.max_virtualspaces(gs))
pspace = Vect[I]((0,0,-P) => 1, (0,0,2*Q-P) => 1, (1,1//2,Q-P) => 1)
n = TensorMap(zeros, elt, pspace ← pspace)
block(n, I((0,0,2*Q-P))) .= 2
block(n, I((1,1//2,Q-P))) .= 1
id = isomorphism(n)
num = sum([correlator(gs, MPSKit.add_util_leg(n), MPSKit.add_util_leg(id), 1, 2),
        correlator(gs, MPSKit.add_util_leg(n), MPSKit.add_util_leg(id), 2, 3),
        correlator(gs, MPSKit.add_util_leg(n), MPSKit.add_util_leg(id), 3, 4),
        correlator(gs, MPSKit.add_util_leg(id), MPSKit.add_util_leg(n), 1, 4)])
println("E0: ", real(E0),"  Particle number: ", real(num))

the results:

[ Info: DMRG2   1:      obj = +5.791932749013e+00       err = 8.7722753525e-01  time = 0.02 sec
[ Info: DMRG2 conv 2:   obj = +5.791932749013e+00       err = 1.3322676296e-15  time = 0.03 sec
E0: 5.7919327490129024  Particle number: 5.0000000000000275
[ Info: DMRG2   1:      obj = +1.353907857862e+01       err = 5.2638077148e-01  time = 0.01 sec
[ Info: DMRG2 conv 2:   obj = +1.353907857862e+01       err = 1.5543122345e-15  time = 0.01 sec
E0: 13.539078578622963  Particle number: 6.00000000000003

[lkdvos]

This is definitely an interesting question, it might actually be worth it to make this a separate issue just to keep it easily accessible in the future.

For the reply, I think the issue is right=Vect[I]((1, 1//2, 1)=> 2) is not entirely what you are thinking: this is an auxiliary space that determines the charge, adding degeneracies to it does not actually alter the charge, you could think about this as now searching for a degenerate subspace instead of a eigenvector. (although the algorithms will probably fail)

What you probably want is right = Vect[I]((0, 0, 2) => 1)) or right = Vect[I]((0, 1, 2) => 1) I think, where this latter one would be the sector of two electrons paired as a triplet. In other words, you want to use the charge of the auxiliary leg, not use the degeneracy.

I haven't tried this out though, so I might be missing something

[ZongYongyue]

This is definitely an interesting question, it might actually be worth it to make this a separate issue just to keep it easily accessible in the future.

For the reply, I think the issue is right=Vect[I]((1, 1//2, 1)=> 2) is not entirely what you are thinking: this is an auxiliary space that determines the charge, adding degeneracies to it does not actually alter the charge, you could think about this as now searching for a degenerate subspace instead of a eigenvector. (although the algorithms will probably fail)

What you probably want is right = Vect[I]((0, 0, 2) => 1)) or right = Vect[I]((0, 1, 2) => 1) I think, where this latter one would be the sector of two electrons paired as a triplet. In other words, you want to use the charge of the auxiliary leg, not use the degeneracy.

I haven't tried this out though, so I might be missing something

Thanks for your correction! I was wrong about charge and degeneracy. If one want add extra particles，it should write as this:

right = Vect[I]([(i, j, k)=>1 for i in ..., for j in ..., for k in ...)

is it right now?

[lkdvos]

I don't think you can combine more than one charge and degeneracy together. What I mean is that you can have a single right = Vect[I]((i, j, k) => 1) and optimize for that, and have to do a separate optimization for a separate sector.

The picture I usually have in my mind is the following: imagining a operator H, this can be fully diagonalized as H * V = V * D with V the square matrix of eigenvectors.
Then, a single eigenvector is obtained by slicing a single column out of V, so this would be a N x 1 matrix. (I'm explicitly keeping the 1-dimension here because of what comes next).

If H was symmetric, it is actually block-diagonal with the blocks labeled by the different sectors.
Similarly, V becomes block diagonal, and you can label the eigenvectors with a sector, so the N x 1 eigenvector actually carries a sector on the 1 as well.
This is precisely the auxiliary space we are adding to the MPS, which is used to select which sector we want to target.
It has to have a single charge and degeneracy 1, since otherwise it would no longer be an eigenvector, but rather an eigenspace that you are targeting (similar to selecting multiple columns from V).
This is not something that our algorithms are set up to handle, and this degeneracy would very likely break DMRG and VUMPS etc, since effectively you would be looking for multiple solutions at the same time, while the algorithms typically assume a unique solution.

This is at least how I try to internalize this, I hope this helps!

[ZongYongyue]

I see——thank you very much for your clear and easy-to-understand explanation!