ensure expval returns normalized values

src/algorithms/expval.jl

@@ -55,38 +55,38 @@ function expectation_value(ψ::AbstractMPS, (inds, O)::Pair)
     # some special cases that avoid using transfer matrices
     if length(sites) == 1
         AC = ψ.AC[sites[1]]
-        return @plansor conj(AC[4 5; 6]) *
-                        conj(Ut[1]) * local_mpo[1][1 5; 3 2] * Ut[2] *
-                        AC[4 3; 6]
-    end
-    if length(sites) == 2 && (sites[1] + 1 == sites[2])
+        E = @plansor conj(AC[4 5; 6]) *
+                     conj(Ut[1]) * local_mpo[1][1 5; 3 2] * Ut[2] *
+                     AC[4 3; 6]
+    elseif length(sites) == 2 && (sites[1] + 1 == sites[2])
         AC = ψ.AC[sites[1]]
         AR = ψ.AR[sites[2]]
         O1, O2 = local_mpo
-        return @plansor conj(AC[4 5; 10]) * conj(Ut[1]) * O1[1 5; 3 8] * AC[4 3; 6] *
-                        conj(AR[10 9; 11]) * Ut[2] * O2[8 9; 7 2] * AR[6 7; 11]
-    end
-
-    # generic case: write as Vl * T^N * Vr
-    # left side
-    T = storagetype(site_type(ψ))
-    @plansor Vl[-1 -2; -3] := isomorphism(T,
-                                          left_virtualspace(ψ, sites[1] - 1),
-                                          left_virtualspace(ψ, sites[1] - 1))[-1; -3] *
-                              conj(Ut[-2])
-
-    # middle
-    M = foldl(zip(sites, local_mpo); init=Vl) do v, (site, o)
-        if o isa Number
-            return scale!(v * TransferMatrix(ψ.AL[site], ψ.AL[site]), o)
-        else
-            return v * TransferMatrix(ψ.AL[site], o, ψ.AL[site])
+        E = @plansor conj(AC[4 5; 10]) * conj(Ut[1]) * O1[1 5; 3 8] * AC[4 3; 6] *
+                     conj(AR[10 9; 11]) * Ut[2] * O2[8 9; 7 2] * AR[6 7; 11]
+    else
+        # generic case: write as Vl * T^N * Vr
+        # left side
+        T = storagetype(site_type(ψ))
+        @plansor Vl[-1 -2; -3] := isomorphism(T,
+                                              left_virtualspace(ψ, sites[1] - 1),
+                                              left_virtualspace(ψ, sites[1] - 1))[-1; -3] *
+                                  conj(Ut[-2])
+
+        # middle
+        M = foldl(zip(sites, local_mpo); init=Vl) do v, (site, o)
+            if o isa Number
+                return scale!(v * TransferMatrix(ψ.AL[site], ψ.AL[site]), o)
+            else
+                return v * TransferMatrix(ψ.AL[site], o, ψ.AL[site])
+            end
         end
+
+        # right side
+        E = @plansor M[1 2; 3] * Ut[2] * ψ.CR[sites[end]][3; 4] *
+                     conj(ψ.CR[sites[end]][1; 4])
     end
 
-    # right side
-    E = @plansor M[1 2; 3] * Ut[2] * ψ.CR[sites[end]][3; 4] *
-                 conj(ψ.CR[sites[end]][1; 4])
     return E / norm(ψ)^2
 end