Hi, ITensor Teams！
For a study of Open Quantum Systems, I need to build L^{\dagger}L where L is an MPO representation of the Liouvillian Superoperator. In the past, I realized that if I simply used LdagL = apply(dag(L),L), I got the wrong result because I was not contracting the output legs of L with the input legs of L^{\dagger}. Due to this, I developed the following function:
using ITensors
using ITensorMPS
function Build_LdagL_MPO(sites, L_MPO::MPO)
Ldag_MPO = dag(L_MPO) #dag results in Index s[j]' (input) and s[j] (output).
for j = 1:length(sites)
Ldag_MPO[j] = Ldag_MPO[j]*delta(sites[j],sites[j]'') #output
Ldag_MPO[j] = Ldag_MPO[j]*delta(sites[j]',sites[j]) #input, will be contracted with L_MPO s[j] if we use apply function.
end
LdagL_MPO = apply(Ldag_MPO, L_MPO) #Index s[j] (input) s[j]''' (output). Let's rename s[j]''' as s[j]' because that is the standard output of an MPO.
for j = 1:length(sites)
LdagL_MPO[j] = LdagL_MPO[j]*delta(sites[j]''',sites[j]')
end
return LdagL_MPO
end
This works well, and I got the correct operator. My problem appears when I try to include symmetries. My L conserves the total number of electrons (not by spin; the following code is a simplified version without dissipation). I wanted to take advantage of that, so I was trying to use:
sites = siteinds("Electron",6, conserve_qns = true)
Physical_site(i::Int64) = 2*i-1
Ancilla_site(i::Int64) = Physical_site(i) + 1
function Build_L_MPO(sites, J::Float64, U::Float64, V::Float64, γ::Float64) #This is a simplified version with only coherent terms
M = length(sites) #Physical + ancilla sites
N = M÷2 #ancilla sites
os = OpSum()
################################################## +H ##################################################
for j=1:N-1
os += -J,"Cdagup",Physical_site(j),"Cup",Physical_site(j+1)
os += -J,"Cdagup",Physical_site(j+1),"Cup",Physical_site(j)
os += -J,"Cdagdn",Physical_site(j),"Cdn",Physical_site(j+1)
os += -J,"Cdagdn",Physical_site(j+1),"Cdn",Physical_site(j)
os += V, "Cdagup * Cup", Physical_site(j), "Cdagup * Cup", Physical_site(j+1)
os += V, "Cdagup * Cup", Physical_site(j), "Cdagdn * Cdn", Physical_site(j+1)
os += V, "Cdagdn * Cdn", Physical_site(j), "Cdagup * Cup", Physical_site(j+1)
os += V, "Cdagdn * Cdn", Physical_site(j), "Cdagdn * Cdn", Physical_site(j+1)
end
for j = 1:N
os += U, "Cdagup * Cup * Cdagdn * Cdn", Physical_site(j)
end
################################################## -H_tilded ##################################################
for j=1:N-1
os -= -J,"Cdagup",Ancilla_site(j),"Cup",Ancilla_site(j+1)
os -= -J,"Cdagup",Ancilla_site(j+1),"Cup",Ancilla_site(j)
os -= -J,"Cdagdn",Ancilla_site(j),"Cdn",Ancilla_site(j+1)
os -= -J,"Cdagdn",Ancilla_site(j+1),"Cdn",Ancilla_site(j)
os -= V, "Cdagup * Cup", Ancilla_site(j), "Cdagup * Cup", Ancilla_site(j+1)
os -= V, "Cdagup * Cup", Ancilla_site(j), "Cdagdn * Cdn", Ancilla_site(j+1)
os -= V, "Cdagdn * Cdn", Ancilla_site(j), "Cdagup * Cup", Ancilla_site(j+1)
os -= V, "Cdagdn * Cdn", Ancilla_site(j), "Cdagdn * Cdn", Ancilla_site(j+1)
end
for j = 1:N
os -= U, "Cdagup * Cup * Cdagdn * Cdn", Ancilla_site(j)
end
L = MPO(os, sites)
return L
end
L_MPO = Build_L_MPO(sites, J, U, V, 0.0); #Everything runs well until here.
LdagL_MPO = Build_LdagL_MPO(sites, L_MPO::MPO)
After the last line I got the following error:
Attempting to contract IndexSet:
((dim=20|id=674|"Link,l=1") <In>
1: QN() => 1
2: QN(("Nf",0,-1),("Sz",0)) => 4
3: QN(("Nf",1,-1),("Sz",-1)) => 2
4: QN(("Nf",-2,-1),("Sz",0)) => 2
5: QN(("Nf",0,-1),("Sz",-2)) => 1
6: QN(("Nf",1,-1),("Sz",1)) => 2
7: QN(("Nf",2,-1),("Sz",0)) => 2
8: QN(("Nf",-1,-1),("Sz",1)) => 2
9: QN(("Nf",0,-1),("Sz",2)) => 1
10: QN(("Nf",-1,-1),("Sz",-1)) => 2
11: QN(("Nf",0,-1),("Sz",0)) => 1, (dim=4|id=95|"Electron,Site,n=1")' <In>
1: QN(("Nf",0,-1),("Sz",0)) => 1
2: QN(("Nf",1,-1),("Sz",1)) => 1
3: QN(("Nf",1,-1),("Sz",-1)) => 1
4: QN(("Nf",2,-1),("Sz",0)) => 1, (dim=4|id=95|"Electron,Site,n=1") <Out>
1: QN(("Nf",0,-1),("Sz",0)) => 1
2: QN(("Nf",1,-1),("Sz",1)) => 1
3: QN(("Nf",1,-1),("Sz",-1)) => 1
4: QN(("Nf",2,-1),("Sz",0)) => 1)
with IndexSet:
((dim=4|id=95|"Electron,Site,n=1") <Out>
1: QN(("Nf",0,-1),("Sz",0)) => 1
2: QN(("Nf",1,-1),("Sz",1)) => 1
3: QN(("Nf",1,-1),("Sz",-1)) => 1
4: QN(("Nf",2,-1),("Sz",0)) => 1, (dim=4|id=95|"Electron,Site,n=1")'' <Out>
1: QN(("Nf",0,-1),("Sz",0)) => 1
2: QN(("Nf",1,-1),("Sz",1)) => 1
3: QN(("Nf",1,-1),("Sz",-1)) => 1
4: QN(("Nf",2,-1),("Sz",0)) => 1)
QN indices must have opposite direction to contract, but indices:
(dim=4|id=95|"Electron,Site,n=1") <Out>
1: QN(("Nf",0,-1),("Sz",0)) => 1
2: QN(("Nf",1,-1),("Sz",1)) => 1
3: QN(("Nf",1,-1),("Sz",-1)) => 1
4: QN(("Nf",2,-1),("Sz",0)) => 1
and:
(dim=4|id=95|"Electron,Site,n=1") <Out>
1: QN(("Nf",0,-1),("Sz",0)) => 1
2: QN(("Nf",1,-1),("Sz",1)) => 1
3: QN(("Nf",1,-1),("Sz",-1)) => 1
4: QN(("Nf",2,-1),("Sz",0)) => 1
do not have opposite directions.
That happens during the first part of my function Build_LdagL_MPO(sites, L_MPO::MPO). So clearly, the function that I created to build L^{\dagger}L is not compatible with an MPO that has Block Structure. Any idea of what I can do? I need L^{\dagger}L as an MPO to run DMRG and find its groundstate (who is equivalent to the eigenstate with 0 eigenvalue of L, i.e., the nonequilibrium steady state). Note that I can not run DMRG directly with L because it is not hermitian and hence its spectrum is complex, and the eigenvalue = 0 is not the minima.
An alternative to find the steady state that I was thinking of is based on the idea that all eigenvalues of L satisfy Re [\lambda_{i}] \leq 0, so the steady state correspond to the eigenvalue with biggest real part. Is there a way to use DMRG with non-Hermitian MPOs and take only the real part of the energy during the optimization? If so, I could use DMRG directly with -L.