Skip to content

Commit 395a3ce

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

File tree

1 file changed

+73
-44
lines changed

1 file changed

+73
-44
lines changed

‎src/algorithms/timestep/wii.jl‎

Lines changed: 73 additions & 44 deletions
Original file line numberDiff line numberDiff line change
@@ -19,62 +19,62 @@ $(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)
3133
init_1 = isometry(storagetype(WD[i]), codomain(WD[i]), domain(WD[i]))
32-
init = [init_1,
34+
init = (init_1,
3335
zero(H[i][1, 1, 1, k]),
3436
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
37+
zero(H[i][j, 1, 1, k]))
38+
39+
y, convhist = exponentiate(1.0, init, exp_alg) do (x1, x2, x3, x4)
40+
@plansor y1[-1 -2; -3 -4] := δ * x1[-1 1; -3 -4] *
41+
H[i][1, 1, 1, end][2 3; 1 4] *
42+
τ[-2 4; 2 3]
43+
44+
@plansor y2[-1 -2; -3 -4] := δ * x2[-1 1; -3 -4] *
45+
H[i][1, 1, 1, end][2 3; 1 4] *
46+
τ[-2 4; 2 3] +
47+
sqrt(δ) *
48+
x1[1 2; -3 4] *
49+
H[i][1, 1, 1, k][-1 -2; 3 -4] *
50+
τ[3 4; 1 2]
51+
52+
@plansor y3[-1 -2; -3 -4] := δ * x3[-1 1; -3 -4] *
53+
H[i][1, 1, 1, end][2 3; 1 4] *
54+
τ[-2 4; 2 3] +
55+
sqrt(δ) *
56+
x1[1 2; -3 4] *
57+
H[i][j, 1, 1, end][-1 -2; 3 -4] *
58+
τ[3 4; 1 2]
59+
60+
@plansor y4[-1 -2; -3 -4] := (δ * x4[-1 1; -3 -4] *
61+
H[i][1, 1, 1, end][2 3; 1 4] *
62+
τ[-2 4; 2 3] +
63+
x1[1 2; -3 4] *
64+
H[i][j, 1, 1, k][-1 -2; 3 -4] *
65+
τ[3 4; 1 2]) +
66+
(sqrt(δ) *
67+
x2[1 2; -3 -4] *
68+
H[i][j, 1, 1, end][-1 -2; 3 4] *
69+
τ[3 4; 1 2]) +
70+
(sqrt(δ) *
71+
x3[-1 4; -3 3] *
72+
H[i][1, 1, 1, k][2 -2; 1 -4] *
73+
τ[3 4; 1 2])
74+
75+
return y1, y2, y3, y4
7776
end
77+
7878
convhist.converged == 0 &&
7979
@warn "failed to exponentiate $(convhist.normres)"
8080

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

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

0 commit comments

Comments
 (0)