Skip to content

Commit bacd649

Browse files
lkdvoslkdvos
committed
WII for FiniteMPOHamiltonian 
1 parent 80cca0f commit bacd649

File tree

1 file changed

+73
-45
lines changed

1 file changed

+73
-45
lines changed

‎src/algorithms/timestep/wii.jl‎

Lines changed: 73 additions & 45 deletions
Original file line numberDiff line numberDiff line change
@@ -19,62 +19,61 @@ $(TYPEDFIELDS)
1919
maxiter::Int = Defaults.maxiter
2020
end
2121

22-
function make_time_mpo(H::InfiniteMPOHamiltonian{T}, dt, alg::WII) where {T}
22+
function make_time_mpo(H::InfiniteMPOHamiltonian, dt::Number, alg::WII)
2323
WA = H.A
2424
WB = H.B
2525
WC = H.C
2626
WD = H.D
2727

2828
δ = dt * (-1im)
29+
exp_alg = Arnoldi(; alg.tol, alg.maxiter)
30+
2931
Wnew = map(1:length(H)) do i
3032
for j in 2:(size(H[i], 1) - 1), k in 2:(size(H[i], 4) - 1)
31-
init_1 = isometry(storagetype(WD[i]), codomain(WD[i]), domain(WD[i]))
32-
init = [init_1,
33+
init = (isometry(storagetype(WD[i]), codomain(WD[i]), domain(WD[i])),
3334
zero(H[i][1, 1, 1, k]),
3435
zero(H[i][j, 1, 1, end]),
35-
zero(H[i][j, 1, 1, k])]
36-
37-
y, convhist = exponentiate(1.0, init,
38-
Arnoldi(; tol=alg.tol, maxiter=alg.maxiter)) do x
39-
out = similar(x)
40-
41-
@plansor out[1][-1 -2; -3 -4] := δ * x[1][-1 1; -3 -4] *
42-
H[i][1, 1, 1, end][2 3; 1 4] *
43-
τ[-2 4; 2 3]
44-
45-
@plansor out[2][-1 -2; -3 -4] := δ * x[2][-1 1; -3 -4] *
46-
H[i][1, 1, 1, end][2 3; 1 4] *
47-
τ[-2 4; 2 3]
48-
@plansor out[2][-1 -2; -3 -4] += sqrt(δ) *
49-
x[1][1 2; -3 4] *
50-
H[i][1, 1, 1, k][-1 -2; 3 -4] *
51-
τ[3 4; 1 2]
52-
53-
@plansor out[3][-1 -2; -3 -4] := δ * x[3][-1 1; -3 -4] *
54-
H[i][1, 1, 1, end][2 3; 1 4] *
55-
τ[-2 4; 2 3]
56-
@plansor out[3][-1 -2; -3 -4] += sqrt(δ) *
57-
x[1][1 2; -3 4] *
58-
H[i][j, 1, 1, end][-1 -2; 3 -4] *
59-
τ[3 4; 1 2]
60-
61-
@plansor out[4][-1 -2; -3 -4] := δ * x[4][-1 1; -3 -4] *
62-
H[i][1, 1, 1, end][2 3; 1 4] *
63-
τ[-2 4; 2 3]
64-
@plansor out[4][-1 -2; -3 -4] += x[1][1 2; -3 4] *
65-
H[i][j, 1, 1, k][-1 -2; 3 -4] *
66-
τ[3 4; 1 2]
67-
@plansor out[4][-1 -2; -3 -4] += sqrt(δ) *
68-
x[2][1 2; -3 -4] *
69-
H[i][j, 1, 1, end][-1 -2; 3 4] *
70-
τ[3 4; 1 2]
71-
@plansor out[4][-1 -2; -3 -4] += sqrt(δ) *
72-
x[3][-1 4; -3 3] *
73-
H[i][1, 1, 1, k][2 -2; 1 -4] *
74-
τ[3 4; 1 2]
75-
76-
return out
36+
zero(H[i][j, 1, 1, k]))
37+
38+
y, convhist = exponentiate(1.0, init, exp_alg) do (x1, x2, x3, x4)
39+
@plansor y1[-1 -2; -3 -4] := δ * x1[-1 1; -3 -4] *
40+
H[i][1, 1, 1, end][2 3; 1 4] *
41+
τ[-2 4; 2 3]
42+
43+
@plansor y2[-1 -2; -3 -4] := δ * x2[-1 1; -3 -4] *
44+
H[i][1, 1, 1, end][2 3; 1 4] *
45+
τ[-2 4; 2 3] +
46+
sqrt(δ) *
47+
x1[1 2; -3 4] *
48+
H[i][1, 1, 1, k][-1 -2; 3 -4] *
49+
τ[3 4; 1 2]
50+
51+
@plansor y3[-1 -2; -3 -4] := δ * x3[-1 1; -3 -4] *
52+
H[i][1, 1, 1, end][2 3; 1 4] *
53+
τ[-2 4; 2 3] +
54+
sqrt(δ) *
55+
x1[1 2; -3 4] *
56+
H[i][j, 1, 1, end][-1 -2; 3 -4] *
57+
τ[3 4; 1 2]
58+
59+
@plansor y4[-1 -2; -3 -4] := (δ * x4[-1 1; -3 -4] *
60+
H[i][1, 1, 1, end][2 3; 1 4] *
61+
τ[-2 4; 2 3] +
62+
x1[1 2; -3 4] *
63+
H[i][j, 1, 1, k][-1 -2; 3 -4] *
64+
τ[3 4; 1 2]) +
65+
(sqrt(δ) *
66+
x2[1 2; -3 -4] *
67+
H[i][j, 1, 1, end][-1 -2; 3 4] *
68+
τ[3 4; 1 2]) +
69+
(sqrt(δ) *
70+
x3[-1 4; -3 3] *
71+
H[i][1, 1, 1, k][2 -2; 1 -4] *
72+
τ[3 4; 1 2])
73+
74+
return y1, y2, y3, y4
7775
end
76+
7877
convhist.converged == 0 &&
7978
@warn "failed to exponentiate $(convhist.normres)"
8079

@@ -99,3 +98,32 @@ function make_time_mpo(H::InfiniteMPOHamiltonian{T}, dt, alg::WII) where {T}
9998

10099
return InfiniteMPO(PeriodicArray(Wnew))
101100
end
101+
102+
# Hack to treat FiniteMPOhamiltonians as Infinite
103+
function make_time_mpo(H::FiniteMPOHamiltonian, dt::Number, alg::WII)
104+
H′ = copy(parent(H))
105+
106+
V_left = left_virtualspace(H[1])
107+
V_left′ = BlockTensorKit.oplus(V_left, oneunit(V_left), oneunit(V_left))
108+
H′[1] = similar(H[1], V_left′ ⊗ space(H[1], 2) ← domain(H[1]))
109+
for (I, v) in nonzero_pairs(H[1])
110+
H′[1][I] = v
111+
end
112+
113+
V_right = right_virtualspace(H[end])
114+
V_right′ = BlockTensorKit.oplus(oneunit(V_right), oneunit(V_right), V_right)
115+
H′[end] = similar(H[end], codomain(H[end]) ← space(H[end], 3)' ⊗ V_right′)
116+
for (I, v) in nonzero_pairs(H[end])
117+
H′[end][I[1], 1, 1, end] = v
118+
end
119+
120+
H′[1][end, 1, 1, end] = H′[1][1, 1, 1, 1]
121+
H′[end][1, 1, 1, 1] = H′[end][end, 1, 1, end]
122+
123+
mpo = make_time_mpo(InfiniteMPOHamiltonian(H′), dt, alg; tol)
124+
125+
# Impose boundary conditions
126+
mpo[1] = mpo[1][1, :, :, :]
127+
mpo[end] = mpo[end][:, :, :, 1]
128+
return remove_orphans(mpo; tol)
129+
end

0 commit comments

Comments
 (0)