|
| 1 | +struct Hᶜ{T} |
| 2 | + ∑h::InfiniteSum{T} |
| 3 | + Hᴸ::InfiniteMPS |
| 4 | + Hᴿ::InfiniteMPS |
| 5 | + ψ::InfiniteCanonicalMPS |
| 6 | + n::Int |
| 7 | +end |
| 8 | + |
| 9 | +struct Hᴬᶜ{T} |
| 10 | + ∑h::InfiniteSum{T} |
| 11 | + Hᴸ::InfiniteMPS |
| 12 | + Hᴿ::InfiniteMPS |
| 13 | + ψ::InfiniteCanonicalMPS |
| 14 | + n::Int |
| 15 | +end |
| 16 | + |
| 17 | +function (H::Hᶜ)(v) |
| 18 | + Hᶜv = H * v |
| 19 | + ## return Hᶜv * δˡ * δʳ |
| 20 | + return noprime(Hᶜv) |
| 21 | +end |
| 22 | + |
| 23 | +function (H::Hᴬᶜ)(v) |
| 24 | + Hᴬᶜv = H * v |
| 25 | + ## return Hᶜv * δˡ⁻¹ * δˢ * δʳ |
| 26 | + return noprime(Hᴬᶜv) |
| 27 | +end |
| 28 | + |
| 29 | +# Struct for use in linear system solver |
| 30 | +struct Aᴸ |
| 31 | + ψ::InfiniteCanonicalMPS |
| 32 | + n::Int |
| 33 | +end |
| 34 | + |
| 35 | +function (A::Aᴸ)(x) |
| 36 | + ψ = A.ψ |
| 37 | + ψᴴ = dag(ψ) |
| 38 | + ψ′ = ψᴴ' |
| 39 | + ψ̃ = prime(linkinds, ψᴴ) |
| 40 | + n = A.n |
| 41 | + |
| 42 | + N = length(ψ) |
| 43 | + #@assert n == N |
| 44 | + |
| 45 | + l = linkinds(only, ψ.AL) |
| 46 | + l′ = linkinds(only, ψ′.AL) |
| 47 | + r = linkinds(only, ψ.AR) |
| 48 | + r′ = linkinds(only, ψ′.AR) |
| 49 | + |
| 50 | + xT = translatecell(x, -1) |
| 51 | + for k in (n - N + 1):n |
| 52 | + xT = xT * ψ.AL[k] * ψ̃.AL[k] |
| 53 | + end |
| 54 | + δˡ = δ(l[n], l′[n]) |
| 55 | + δʳ = δ(r[n], r′[n]) |
| 56 | + xR = x * ψ.C[n] * ψ′.C[n] * dag(δʳ) * denseblocks(δˡ) |
| 57 | + return xT - xR |
| 58 | +end |
| 59 | + |
| 60 | +function left_environment(hᴸ, 𝕙ᴸ, ψ; tol=1e-15) |
| 61 | + ψ̃ = prime(linkinds, dag(ψ)) |
| 62 | + N = nsites(ψ) |
| 63 | + |
| 64 | + Aᴺ = Aᴸ(ψ, N) |
| 65 | + Hᴸᴺ¹, info = linsolve(Aᴺ, 𝕙ᴸ[N], 1, -1; tol=tol) |
| 66 | + # Get the rest of the environments in the unit cell |
| 67 | + Hᴸ = InfiniteMPS(Vector{ITensor}(undef, N)) |
| 68 | + Hᴸ[N] = Hᴸᴺ¹ |
| 69 | + for n in 1:(N - 1) |
| 70 | + Hᴸ[n] = Hᴸ[n - 1] * ψ.AL[n] * ψ̃.AL[n] + hᴸ[n] |
| 71 | + end |
| 72 | + return Hᴸ |
| 73 | +end |
| 74 | + |
| 75 | +# Struct for use in linear system solver |
| 76 | +struct Aᴿ |
| 77 | + hᴿ::InfiniteMPS |
| 78 | + ψ::InfiniteCanonicalMPS |
| 79 | + n::Int |
| 80 | +end |
| 81 | + |
| 82 | +function (A::Aᴿ)(x) |
| 83 | + hᴿ = A.hᴿ |
| 84 | + ψ = A.ψ |
| 85 | + ψᴴ = dag(ψ) |
| 86 | + ψ′ = ψᴴ' |
| 87 | + ψ̃ = prime(linkinds, ψᴴ) |
| 88 | + n = A.n |
| 89 | + |
| 90 | + N = length(ψ) |
| 91 | + @assert n == N |
| 92 | + |
| 93 | + l = linkinds(only, ψ.AL) |
| 94 | + l′ = linkinds(only, ψ′.AL) |
| 95 | + r = linkinds(only, ψ.AR) |
| 96 | + r′ = linkinds(only, ψ′.AR) |
| 97 | + |
| 98 | + xT = x |
| 99 | + for k in reverse(1:N) |
| 100 | + xT = xT * ψ.AR[k] * ψ̃.AR[k] |
| 101 | + end |
| 102 | + xT = translatecell(xT, 1) |
| 103 | + δˡ = δ(l[n], l′[n]) |
| 104 | + δʳ = δ(r[n], r′[n]) |
| 105 | + xR = x * ψ.C[n] * ψ′.C[n] * δˡ * denseblocks(dag(δʳ)) |
| 106 | + return xT - xR |
| 107 | +end |
| 108 | + |
| 109 | +function right_environment(hᴿ, ψ; tol=1e-15) |
| 110 | + ψ̃ = prime(linkinds, dag(ψ)) |
| 111 | + # XXX: replace with `nsites` |
| 112 | + #N = nsites(ψ) |
| 113 | + N = length(ψ) |
| 114 | + |
| 115 | + A = Aᴿ(hᴿ, ψ, N) |
| 116 | + Hᴿᴺ¹, info = linsolve(A, hᴿ[N], 1, -1; tol=tol) |
| 117 | + # Get the rest of the environments in the unit cell |
| 118 | + Hᴿ = InfiniteMPS(Vector{ITensor}(undef, N)) |
| 119 | + Hᴿ[N] = Hᴿᴺ¹ |
| 120 | + for n in reverse(1:(N - 1)) |
| 121 | + Hᴿ[n] = Hᴿ[n + 1] * ψ.AR[n + 1] * ψ̃.AR[n + 1] + hᴿ[n] |
| 122 | + end |
| 123 | + return Hᴿ |
| 124 | +end |
| 125 | + |
| 126 | +# TODO Generate all environments, why? Only one is needed in the sequential version |
| 127 | +function right_environment(hᴿ, 𝕙ᴿ, ψ; tol=1e-15) |
| 128 | + ψ̃ = prime(linkinds, dag(ψ)) |
| 129 | + N = nsites(ψ) |
| 130 | + |
| 131 | + A = Aᴿ(hᴿ, ψ, N) |
| 132 | + Hᴿᴺ¹, info = linsolve(A, 𝕙ᴿ[N], 1, -1; tol=tol) |
| 133 | + # Get the rest of the environments in the unit cell |
| 134 | + Hᴿ = InfiniteMPS(Vector{ITensor}(undef, N)) |
| 135 | + Hᴿ[N] = Hᴿᴺ¹ |
| 136 | + for n in reverse(1:(N - 1)) |
| 137 | + Hᴿ[n] = Hᴿ[n + 1] * ψ.AR[n + 1] * ψ̃.AR[n + 1] + hᴿ[n] |
| 138 | + end |
| 139 | + return Hᴿ |
| 140 | +end |
| 141 | + |
| 142 | +function tdvp_iteration(args...; multisite_update_alg="sequential", kwargs...) |
| 143 | + if multisite_update_alg == "sequential" |
| 144 | + return tdvp_iteration_sequential(args...; kwargs...) |
| 145 | + elseif multisite_update_alg == "parallel" |
| 146 | + return tdvp_iteration_parallel(args...; kwargs...) |
| 147 | + else |
| 148 | + error( |
| 149 | + "Multisite update algorithm multisite_update_alg = $multisite_update_alg not supported, use \"parallel\" or \"sequential\"", |
| 150 | + ) |
| 151 | + end |
| 152 | +end |
| 153 | + |
| 154 | +function tdvp( |
| 155 | + solver::Function, |
| 156 | + ∑h, |
| 157 | + ψ; |
| 158 | + maxiter=10, |
| 159 | + tol=1e-8, |
| 160 | + outputlevel=1, |
| 161 | + multisite_update_alg="sequential", |
| 162 | + time_step, |
| 163 | + solver_tol=(x -> x / 100), |
| 164 | +) |
| 165 | + N = nsites(ψ) |
| 166 | + (ϵᴸ!) = fill(tol, nsites(ψ)) |
| 167 | + (ϵᴿ!) = fill(tol, nsites(ψ)) |
| 168 | + if outputlevel > 0 |
| 169 | + println("Running VUMPS with multisite_update_alg = $multisite_update_alg") |
| 170 | + flush(stdout) |
| 171 | + flush(stderr) |
| 172 | + end |
| 173 | + for iter in 1:maxiter |
| 174 | + ψ, (eᴸ, eᴿ) = tdvp_iteration( |
| 175 | + solver, |
| 176 | + ∑h, |
| 177 | + ψ; |
| 178 | + (ϵᴸ!)=(ϵᴸ!), |
| 179 | + (ϵᴿ!)=(ϵᴿ!), |
| 180 | + multisite_update_alg=multisite_update_alg, |
| 181 | + solver_tol=solver_tol, |
| 182 | + time_step=time_step, |
| 183 | + ) |
| 184 | + ϵᵖʳᵉˢ = max(maximum(ϵᴸ!), maximum(ϵᴿ!)) |
| 185 | + maxdimψ = maxlinkdim(ψ[0:(N + 1)]) |
| 186 | + if outputlevel > 0 |
| 187 | + println( |
| 188 | + "VUMPS iteration $iter (out of maximum $maxiter). Bond dimension = $maxdimψ, energy = $((eᴸ, eᴿ)), ϵᵖʳᵉˢ = $ϵᵖʳᵉˢ, tol = $tol", |
| 189 | + ) |
| 190 | + flush(stdout) |
| 191 | + flush(stderr) |
| 192 | + end |
| 193 | + if ϵᵖʳᵉˢ < tol |
| 194 | + println( |
| 195 | + "Precision error $ϵᵖʳᵉˢ reached tolerance $tol, stopping VUMPS after $iter iterations (of a maximum $maxiter).", |
| 196 | + ) |
| 197 | + flush(stdout) |
| 198 | + flush(stderr) |
| 199 | + break |
| 200 | + end |
| 201 | + end |
| 202 | + return ψ |
| 203 | +end |
| 204 | + |
| 205 | +function vumps_solver(M, time_step, v₀, solver_tol) |
| 206 | + λ⃗, v⃗, info = eigsolve(M, v₀, 1, :SR; ishermitian=true, tol=solver_tol) |
| 207 | + return λ⃗[1], v⃗[1], info |
| 208 | +end |
| 209 | + |
| 210 | +return function tdvp_solver(M, time_step, v₀, solver_tol) |
| 211 | + v, info = exponentiate(M, time_step, v₀; ishermitian=true, tol=solver_tol) |
| 212 | + v = v / norm(v) |
| 213 | + return nothing, v, info |
| 214 | +end |
| 215 | + |
| 216 | +function vumps( |
| 217 | + args...; time_step=-Inf, eigsolve_tol=(x -> x / 100), solver_tol=eigsolve_tol, kwargs... |
| 218 | +) |
| 219 | + @assert isinf(time_step) && time_step < 0 |
| 220 | + println("Using VUMPS solver with time step $time_step") |
| 221 | + flush(stdout) |
| 222 | + flush(stderr) |
| 223 | + return tdvp(vumps_solver, args...; time_step=time_step, solver_tol=solver_tol, kwargs...) |
| 224 | +end |
| 225 | + |
| 226 | +function tdvp(args...; time_step, solver_tol=(x -> x / 100), kwargs...) |
| 227 | + solver = if !isinf(time_step) |
| 228 | + println("Using TDVP solver with time step $time_step") |
| 229 | + flush(stdout) |
| 230 | + flush(stderr) |
| 231 | + tdvp_solver |
| 232 | + elseif time_step < 0 |
| 233 | + # Call VUMPS instead |
| 234 | + println("Using VUMPS solver with time step $time_step") |
| 235 | + flush(stdout) |
| 236 | + flush(stderr) |
| 237 | + vumps_solver |
| 238 | + else |
| 239 | + error("Time step $time_step not supported.") |
| 240 | + end |
| 241 | + return tdvp(solver, args...; time_step=time_step, solver_tol=solver_tol, kwargs...) |
| 242 | +end |
0 commit comments