diff --git a/src/algorithms/expval.jl b/src/algorithms/expval.jl index 120f9880..fd0f6f27 100644 --- a/src/algorithms/expval.jl +++ b/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