Quasiparticle type stability improvements

Test type stability excitations

Author: lkdvos

src/algorithms/excitation/quasiparticleexcitation.jl

@@ -342,52 +342,19 @@ function _effective_excitation_local_apply(site, ϕ,
     return B′
 end
 
-function effective_excitation_renormalization_energy(H::Union{InfiniteMPO,
-                                                              InfiniteMPOHamiltonian,
-                                                              FiniteMPOHamiltonian}, ϕ,
-                                                     lenvs,
-                                                     renvs)
+function effective_excitation_renormalization_energy(H, ϕ, lenvs, renvs)
     ψ_left = ϕ.left_gs
-    E_left = map(1:length(ϕ)) do site
-        return contract_mpo_expval(ψ_left.AC[site], leftenv(lenvs, site, ψ_left),
-                                   H[site], rightenv(lenvs, site, ψ_left))
-    end
-
-    ϕ.trivial && return E_left
-
     ψ_right = ϕ.right_gs
-    E_right = map(1:length(ϕ)) do site
-        return contract_mpo_expval(ψ_right.AC[site], leftenv(renvs, site, ψ_right),
-                                   H[site], rightenv(renvs, site, ψ_right))
-    end
-
-    return (E_left .+ E_right) ./ 2
-end
-
-function effective_excitation_renormalization_energy(H, ϕ, lenvs, renvs)
-    E_left = map(1:length(ϕ)) do loc
-        AC = ϕ.left_gs.AC[loc]
-        GL = leftenv(lenvs, loc, ϕ.left_gs)
-        GR = rightenv(lenvs, loc, ϕ.left_gs)
-        return sum(keys(H[loc]); init=zero(scalartype(ϕ))) do (j, k)
-            return @plansor conj(AC[2 6; 4]) * GL[j][2 5; 3] * AC[3 7; 1] *
-                            H[loc][j, k][5 6; 7 8] *
-                            GR[k][1 8; 4]
+    E = Vector{scalartype(ϕ)}(undef, length(ϕ))
+    for i in eachindex(E)
+        E[i] = contract_mpo_expval(ψ_left.AC[i], leftenv(lenvs, i, ψ_left),
+                                   H[i], rightenv(lenvs, i, ψ_left))
+        if !ϕ.trivial
+            E[i] += contract_mpo_expval(ψ_right.AC[i], leftenv(renvs, i, ψ_right),
+                                        H[i], rightenv(renvs, i, ψ_right))
+            E[i] /= 2
         end
     end
-
-    ϕ.trivial && return E_left
-
-    E_right = map(1:length(ϕ)) do loc
-        AC = ϕ.right_gs.AC[loc]
-        GL = leftenv(renvs, loc, ϕ.right_gs)
-        GR = rightenv(renvs, loc, ϕ.right_gs)
-        return sum(keys(H[loc]); init=zero(scalartype(ϕ))) do (j, k)
-            return @plansor conj(AC[2 6; 4]) * GL[j][2 5; 3] * AC[3 7; 1] *
-                            H[loc][j, k][5 6; 7 8] *
-                            GR[k][1 8; 4]
-        end
-    end
-
-    return (E_left .+ E_right) ./ 2
+    return E
 end
+

src/states/quasiparticle_state.jl

@@ -59,13 +59,11 @@ function RightGaugedQP(datfun, left_gs, right_gs=left_gs;
 end
 
 #gauge dependent code
-function Base.similar(v::LeftGaugedQP, T=scalartype(v))
-    return LeftGaugedQP(v.left_gs, v.right_gs, v.VLs, map(e -> similar(e, T), v.Xs),
-                        v.momentum)
+function Base.similar(v::LeftGaugedQP, ::Type{T}=scalartype(v)) where {T<:Number}
+    return LeftGaugedQP(v.left_gs, v.right_gs, v.VLs, similar.(v.Xs, T), v.momentum)
 end
-function Base.similar(v::RightGaugedQP, T=scalartype(v))
-    return RightGaugedQP(v.left_gs, v.right_gs, map(e -> similar(e, T), v.Xs), v.VRs,
-                         v.momentum)
+function Base.similar(v::RightGaugedQP, ::Type{T}=scalartype(v)) where {T<:Number}
+    return RightGaugedQP(v.left_gs, v.right_gs, similar.(v.Xs, T), v.VRs, v.momentum)
 end
 
 Base.getindex(v::LeftGaugedQP, i::Int) = v.VLs[mod1(i, end)] * v.Xs[mod1(i, end)];

test/algorithms.jl

@@ -435,7 +435,8 @@ end
         ψ = InfiniteMPS([ℙ^3, ℙ^3], [ℙ^48, ℙ^48])
         ψ, envs, _ = find_groundstate(ψ, H; maxiter=400, verbosity=verbosity_conv,
                                       tol=1e-10)
-        energies, ϕs = excitations(H, QuasiparticleAnsatz(), Float64(pi), ψ, envs)
+        energies, ϕs = @constinferred excitations(H, QuasiparticleAnsatz(), Float64(pi), ψ,
+                                                  envs)
         @test energies[1] ≈ 0.41047925 atol = 1e-4
         @test variance(ϕs[1], H) < 1e-8
     end
@@ -445,7 +446,8 @@ end
         ψ, envs, _ = leading_boundary(ψ, H,
                                       VUMPS(; maxiter=400, verbosity=verbosity_conv,
                                             tol=1e-10))
-        energies, ϕs = excitations(H, QuasiparticleAnsatz(), [0.0, Float64(pi / 2)], ψ,
+        energies, ϕs = @constinferred excitations(H, QuasiparticleAnsatz(),
+                                                  [0.0, Float64(pi / 2)], ψ,
                                    envs; verbosity=0)
         @test abs(energies[1]) > abs(energies[2]) # has a minimum at pi/2
     end
@@ -464,20 +466,20 @@ end
             ψ, envs, = find_groundstate(ψ, H; verbosity)
 
             # find energy with quasiparticle ansatz
-            energies_QP, ϕs = excitations(H, QuasiparticleAnsatz(), ψ, envs)
+            energies_QP, ϕs = @constinferred excitations(H, QuasiparticleAnsatz(), ψ, envs)
             @test variance(ϕs[1], H) < 1e-6
 
             # find energy with normal dmrg
             for gsalg in (DMRG(; verbosity, tol=1e-6),
                           DMRG2(; verbosity, tol=1e-6, trscheme=truncbelow(1e-4)))
-                energies_dm, _ = excitations(H, FiniteExcited(; gsalg), ψ)
+                energies_dm, _ = @constinferred excitations(H, FiniteExcited(; gsalg), ψ)
                 @test energies_dm[1] ≈ energies_QP[1] + expectation_value(ψ, H, envs) atol = 1e-4
             end
 
             # find energy with Chepiga ansatz
-            energies_ch, _ = excitations(H, ChepigaAnsatz(), ψ, envs)
+            energies_ch, _ = @constinferred excitations(H, ChepigaAnsatz(), ψ, envs)
             @test energies_ch[1] ≈ energies_QP[1] + expectation_value(ψ, H, envs) atol = 1e-4
-            energies_ch2, _ = excitations(H, ChepigaAnsatz2(), ψ, envs)
+            energies_ch2, _ = @constinferred excitations(H, ChepigaAnsatz2(), ψ, envs)
             @test energies_ch2[1] ≈ energies_QP[1] + expectation_value(ψ, H, envs) atol = 1e-4
             return energies_QP[1]
         end