Skip to content

Commit 1f19fd1

Browse files
JanReimersJanReimers
authored
InfiniteMPOMatrix to InfiniteMPO conversion (#76)
1 parent 61d5bbb commit 1f19fd1

File tree

4 files changed

+260
-1
lines changed

4 files changed

+260
-1
lines changed

‎src/infinitempo.jl‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -11,3 +11,5 @@ mutable struct InfiniteMPO <: AbstractInfiniteMPS
1111
end
1212

1313
translator(mpo::InfiniteMPO) = mpo.data.translator
14+
15+
InfiniteMPO(data::CelledVector{ITensor}) = InfiniteMPO(data, 0, size(data, 1), false)

‎src/infinitempomatrix.jl‎

Lines changed: 109 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -7,6 +7,7 @@ mutable struct InfiniteMPOMatrix <: AbstractInfiniteMPS
77
end
88

99
translator(mpo::InfiniteMPOMatrix) = mpo.data.translator
10+
data(mpo::InfiniteMPOMatrix) = mpo.data
1011

1112
# TODO better printing?
1213
function Base.show(io::IO, M::InfiniteMPOMatrix)
@@ -43,4 +44,111 @@ function InfiniteMPOMatrix(data::Vector{Matrix{ITensor}}, translator::Function)
4344
return InfiniteMPOMatrix(CelledVector(data, translator), 0, size(data)[1], false)
4445
end
4546

46-
#nrange(H::InfiniteMPOMatrix) = [size(H[j])[1] - 1 for j in 1:nsites(H)]
47+
#
48+
# Hm should have the form below. Only link indices are shown.
49+
#
50+
# I 0 0 0 0
51+
# M-->l=n 0 0 0 0
52+
# 0 l=n-->M-->l=n-1 ... 0 0 0
53+
# : : : : :
54+
# 0 0 ... l=2-->M-->l=1 0 0
55+
# 0 0 ... 0 l=1-->M I
56+
#
57+
# We need to capture the corner I's and the sub-diagonal Ms, and paste them together into an ITensor.
58+
# This is facilutated by making all elements of Hm into order(4) tensors by adding dummy Dw=1 indices.
59+
# We ITensors.directsum() to join all the blocks into the final ITensor.
60+
# This code should be 100% dense/blocks-sparse agnostic.
61+
#
62+
function cat_to_itensor(Hm::Matrix{ITensor}, inds_array)::Tuple{ITensor,Index,Index}
63+
lx, ly = size(Hm)
64+
@assert lx == ly
65+
#
66+
# Extract the sub diagonal
67+
#
68+
Ms = map(n -> Hm[lx - n + 1, ly - n], 1:(lx - 1)) #Get the sub diagonal into an array.
69+
N = length(Ms)
70+
@assert N == lx - 1
71+
#
72+
# Convert edge Ms to order 4 ITensors using the dummy index.
73+
#
74+
il0 = inds_array[1][1]
75+
T = eltype(Ms[1])
76+
Ms[1] *= onehot(T, il0 => 1)
77+
Ms[N] *= onehot(T, dag(il0) => 1) #We don't need distinct index IDs for this. directsum makes all new index anyway.
78+
#
79+
# Start direct sums to build the single ITensor. Order is critical to get the desired form.
80+
#
81+
H = Hm[1, 1] * onehot(T, il0' => 1) * onehot(T, dag(il0) => 1) #Bootstrap with the I op in the top left of Hm
82+
H, il = directsum(H => il0', Ms[N] => inds_array[N][1]) #1-D directsum to put M[N] directly below the I op
83+
ir = dag(il0) #setup recusion.
84+
for i in (N - 1):-1:1 #2-D directsums to place new blocks below and to the right.
85+
H, (il, ir) = directsum(
86+
H => (il, ir), Ms[i] => inds_array[i]; tags=["Link,left", "Link,right"]
87+
)
88+
end
89+
IN = Hm[N + 1, N + 1] * onehot(T, dag(il0) => 1) * onehot(T, il => dim(il)) #I op in the bottom left of Hm
90+
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]
91+
92+
return H, il, ir
93+
end
94+
95+
function find_all_links(Hms::InfiniteMPOMatrix)
96+
is = inds(Hms[1][1, 1]) #site inds
97+
ir = noncommonind(Hms[1][2, 1], is) #This op should have one link index pointing to the right.
98+
#
99+
# Set up return array of 2-tuples
100+
#
101+
indexT = typeof(ir)
102+
TupleT = NTuple{2,indexT}
103+
inds_array = Vector{TupleT}[]
104+
#
105+
# Make a dummy index
106+
#
107+
il0 = Index(ITensors.trivial_space(ir); dir=dir(dag(ir)), tags="Link,l=0")
108+
#
109+
# Loop over sites and nonzero matrix elements which are linked into the next
110+
# and previous iMPOMatrix.
111+
#
112+
for n in 1:nsites(Hms)
113+
Hm = Hms[n]
114+
inds_n = TupleT[]
115+
lx, ly = size(Hm)
116+
@assert lx == ly
117+
for iy in (ly - 1):-1:1
118+
ix = iy + 1
119+
il = ix < lx ? commonind(Hm[ix, iy], Hms[n - 1][ix + 1, iy + 1]) : dag(il0)
120+
ir = iy > 1 ? commonind(Hm[ix, iy], Hms[n + 1][ix - 1, iy - 1]) : il0
121+
push!(inds_n, (il, ir))
122+
end
123+
push!(inds_array, inds_n)
124+
end
125+
return inds_array
126+
end
127+
128+
#
129+
# Hm is the InfiniteMPOMatrix
130+
# Hlrs is an array of {ITensor,Index,Index}s, one for each site in the unit cell.
131+
# Hi is a CelledVector of ITensors.
132+
#
133+
function InfiniteMPO(Hm::InfiniteMPOMatrix)
134+
inds_array = find_all_links(Hm)
135+
Hlrs = cat_to_itensor.(Hm, inds_array) #return an array of {ITensor,Index,Index}
136+
#
137+
# Unpack the array of tuples into three arrays. And also get an array site indices.
138+
#
139+
Hi = CelledVector([Hlr[1] for Hlr in Hlrs], translator(Hm))
140+
ils = CelledVector([Hlr[2] for Hlr in Hlrs], translator(Hm))
141+
irs = CelledVector([Hlr[3] for Hlr in Hlrs], translator(Hm))
142+
site_inds = [commoninds(Hm[j][1, 1], Hm[j][end, end])[1] for j in 1:nsites(Hm)]
143+
#
144+
# Create new tags with proper cell and link numbers. Also daisy chain
145+
# all the indices so right index at j = dag(left index at j+1)
146+
#
147+
for j in 1:nsites(Hm)
148+
newTag = "Link,c=$(getcell(site_inds[j])),l=$(getsite(site_inds[j]))"
149+
ir = replacetags(irs[j], tags(irs[j]), newTag) #new right index
150+
Hi[j] = replaceinds(Hi[j], [irs[j]], [ir])
151+
Hi[j + 1] = replaceinds(Hi[j + 1], [ils[j + 1]], [dag(ir)])
152+
end
153+
return InfiniteMPO(Hi)
154+
end

‎src/models/models.jl‎

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -82,6 +82,14 @@ end
8282
import ITensors: op
8383
op(::OpName"Zero", ::SiteType, s::Index) = ITensor(s', dag(s))
8484

85+
function InfiniteMPO(model::Model, s::CelledVector; kwargs...)
86+
return InfiniteMPO(model, s, translator(s); kwargs...)
87+
end
88+
89+
function InfiniteMPO(model::Model, s::CelledVector, translator::Function; kwargs...)
90+
return InfiniteMPO(InfiniteMPOMatrix(model, s, translator; kwargs...))
91+
end
92+
8593
function InfiniteMPOMatrix(model::Model, s::CelledVector; kwargs...)
8694
return InfiniteMPOMatrix(model, s, translator(s); kwargs...)
8795
end

‎test/test_iMPOConversions.jl‎

Lines changed: 141 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,141 @@
1+
using ITensors
2+
using ITensorInfiniteMPS
3+
using Test
4+
#
5+
# InfiniteMPO has dangling links at the end of the chain. We contract these on the outside
6+
# with l,r terminating vectors, to make a finite lattice MPO.
7+
#
8+
function terminate(h::InfiniteMPO)::MPO
9+
Ncell = nsites(h)
10+
# left termination vector
11+
il0 = commonind(h[1], h[0])
12+
l = ITensor(0.0, il0)
13+
l[il0 => dim(il0)] = 1.0 #assuming lower reg form in h
14+
# right termination vector
15+
iln = commonind(h[Ncell], h[Ncell + 1])
16+
r = ITensor(0.0, iln)
17+
r[iln => 1] = 1.0 #assuming lower reg form in h
18+
# build up a finite MPO
19+
hf = MPO(Ncell)
20+
hf[1] = dag(l) * h[1] #left terminate
21+
hf[Ncell] = h[Ncell] * dag(r) #right terminate
22+
for n in 2:(Ncell - 1)
23+
hf[n] = h[n] #fill in the bulk.
24+
end
25+
return hf
26+
end
27+
#
28+
# Terminate and then call expect
29+
# for inf ψ and finite h, which is already supported in src/infinitecanonicalmps.jl
30+
#
31+
function ITensors.expect(ψ::InfiniteCanonicalMPS, h::InfiniteMPO)
32+
return expect(ψ, terminate(h)) #defer to src/infinitecanonicalmps.jl
33+
end
34+
35+
#H = ΣⱼΣn (½ S⁺ⱼS⁻ⱼ₊n + ½ S⁻ⱼS⁺ⱼ₊n + SᶻⱼSᶻⱼ₊n)
36+
function ITensorInfiniteMPS.unit_cell_terms(::Model"heisenbergNNN"; NNN::Int64)
37+
opsum = OpSum()
38+
for n in 1:NNN
39+
J = 1.0 / n
40+
opsum += J * 0.5, "S+", 1, "S-", 1 + n
41+
opsum += J * 0.5, "S-", 1, "S+", 1 + n
42+
opsum += J, "Sz", 1, "Sz", 1 + n
43+
end
44+
return opsum
45+
end
46+
47+
function ITensorInfiniteMPS.unit_cell_terms(::Model"hubbardNNN"; NNN::Int64)
48+
U::Float64 = 0.25
49+
t::Float64 = 1.0
50+
V::Float64 = 0.5
51+
opsum = OpSum()
52+
opsum += (U, "Nupdn", 1)
53+
for n in 1:NNN
54+
tj, Vj = t / n, V / n
55+
opsum += -tj, "Cdagup", 1, "Cup", 1 + n
56+
opsum += -tj, "Cdagup", 1 + n, "Cup", 1
57+
opsum += -tj, "Cdagdn", 1, "Cdn", 1 + n
58+
opsum += -tj, "Cdagdn", 1 + n, "Cdn", 1
59+
opsum += Vj, "Ntot", 1, "Ntot", 1 + n
60+
end
61+
return opsum
62+
end
63+
64+
function ITensors.space(::SiteType"FermionK", pos::Int; p=1, q=1, conserve_momentum=true)
65+
if !conserve_momentum
66+
return [QN("Nf", -p) => 1, QN("Nf", q - p) => 1]
67+
else
68+
return [
69+
QN(("Nf", -p), ("NfMom", -p * pos)) => 1,
70+
QN(("Nf", q - p), ("NfMom", (q - p) * pos)) => 1,
71+
]
72+
end
73+
end
74+
75+
function fermion_momentum_translator(i::Index, n::Integer; N=6)
76+
#@show n
77+
ts = tags(i)
78+
translated_ts = ITensorInfiniteMPS.translatecelltags(ts, n)
79+
new_i = replacetags(i, ts => translated_ts)
80+
for j in 1:length(new_i.space)
81+
ch = new_i.space[j][1][1].val
82+
mom = new_i.space[j][1][2].val
83+
new_i.space[j] = Pair(QN(("Nf", ch), ("NfMom", mom + n * N * ch)), new_i.space[j][2])
84+
end
85+
return new_i
86+
end
87+
88+
function ITensors.op!(Op::ITensor, opname::OpName, ::SiteType"FermionK", s::Index...)
89+
return ITensors.op!(Op, opname, SiteType("Fermion"), s...)
90+
end
91+
92+
@testset verbose = true "InfiniteMPOMatrix -> InfiniteMPO" begin
93+
ferro(n) = "↑"
94+
antiferro(n) = isodd(n) ? "↑" : "↓"
95+
96+
models = [(Model"heisenbergNNN"(), "S=1/2"), (Model"hubbardNNN"(), "Electron")]
97+
@testset "H=$model, Ncell=$Ncell, NNN=$NNN, Antiferro=$Af, qns=$qns" for (model, site) in
98+
models,
99+
qns in [false, true],
100+
Ncell in 2:6,
101+
NNN in 1:(Ncell - 1),
102+
Af in [true, false]
103+
104+
if isodd(Ncell) && Af #skip test since Af state does fit inside odd cells.
105+
continue
106+
end
107+
initstate(n) = Af ? antiferro(n) : ferro(n)
108+
model_kwargs = (NNN=NNN,)
109+
s = infsiteinds(site, Ncell; initstate, conserve_qns=qns)
110+
ψ = InfMPS(s, initstate)
111+
Hi = InfiniteMPO(model, s; model_kwargs...)
112+
Hs = InfiniteSum{MPO}(model, s; model_kwargs...)
113+
Es = expect(ψ, Hs)
114+
Ei = expect(ψ, Hi)
115+
#@show Es Ei
116+
@test sum(Es[1:(Ncell - NNN)]) ≈ Ei atol = 1e-14
117+
end
118+
119+
@testset "FQHE Hamitonian" begin
120+
N = 6
121+
model = Model"fqhe_2b_pot"()
122+
model_params = (Vs=[1.0, 0.0, 1.0, 0.0, 0.1], Ly=3.0, prec=1e-8)
123+
trf = fermion_momentum_translator
124+
function initstate(n)
125+
mod1(n, 3) == 1 && return 2
126+
return 1
127+
end
128+
p = 1
129+
q = 3
130+
conserve_momentum = true
131+
s = infsiteinds("FermionK", N; translator=trf, initstate, conserve_momentum, p, q)
132+
ψ = InfMPS(s, initstate)
133+
Hs = InfiniteSum{MPO}(model, s; model_params...)
134+
Hi = InfiniteMPO(model, s, trf; model_params...)
135+
Es = expect(ψ, Hs)
136+
Ei = expect(ψ, Hi)
137+
#@show Es Ei
138+
@test Es[1] ≈ Ei atol = 1e-14
139+
end
140+
end
141+
nothing

0 commit comments

Comments
 (0)