@@ -37,53 +37,42 @@ function expectation_value end
3737function expectation_value (ψ:: AbstractMPS , (inds, O):: Pair )
3838 @boundscheck foreach (Base. Fix1 (Base. checkbounds, ψ), inds)
3939
40- sites, local_mpo = instantiate_operator (physicalspace (ψ), inds => O)
41- @assert _firstspace (first (local_mpo)) == oneunit (_firstspace (first (local_mpo))) ==
42- dual (_lastspace (last (local_mpo)))
43- for (site, o) in zip (sites, local_mpo)
44- if o isa MPOTensor
45- physicalspace (ψ, site) == physicalspace (o) ||
46- throw (SpaceMismatch (" physical space does not match at site $site "
47- end
40+ # special cases that can be handled more efficiently
41+ if length (inds) == 1
42+ return local_expectation_value1 (ψ, inds[1 ], O)
43+ elseif length (inds) == 2 && eltype (inds) <: Integer && inds[1 ] + 1 == inds[2 ]
44+ return local_expectation_value2 (ψ, inds[1 ], O)
4845 end
4946
50- Ut = fill_data! (similar (local_mpo[1 ], _firstspace (first (local_mpo))), one)
51-
52- # some special cases that avoid using transfer matrices
53- if length (sites) == 1
54- AC = ψ. AC[sites[1 ]]
55- E = @plansor conj (AC[4 5 ; 6 ]) * conj (Ut[1 ]) * local_mpo[1 ][1 5 ; 3 2 ] * Ut[2 ] *
56- AC[4 3 ; 6 ]
57- elseif length (sites) == 2 && (sites[1 ] + 1 == sites[2 ])
58- AC = ψ. AC[sites[1 ]]
59- AR = ψ. AR[sites[2 ]]
60- O1, O2 = local_mpo
61- E = @plansor conj (AC[4 5 ; 10 ]) * conj (Ut[1 ]) * O1[1 5 ; 3 8 ] * AC[4 3 ; 6 ] *
62- conj (AR[10 9 ; 11 ]) * Ut[2 ] * O2[8 9 ; 7 2 ] * AR[6 7 ; 11 ]
63- else
64- # generic case: write as Vl * T^N * Vr
65- # left side
66- T = storagetype (site_type (ψ))
67- @plansor Vl[- 1 - 2 ; - 3 ] := isomorphism (
68- T, left_virtualspace (ψ, sites[1 ]), left_virtualspace (ψ, sites[1 ])
69- )[- 1 ; - 3 ] *
70- conj (Ut[- 2 ])
71-
72- # middle
73- M = foldl (zip (sites, local_mpo); init = Vl) do v, (site, o)
74- if o isa Number
75- return scale! (v * TransferMatrix (ψ. AL[site], ψ. AL[site]), o)
76- else
77- return v * TransferMatrix (ψ. AL[site], o, ψ. AL[site])
78- end
79- end
47+ # generic case: convert to MPO and write as Vl * T^N * Vr
48+ sites, local_mpo = instantiate_operator (ψ, inds => O)
49+
50+ # left side
51+ Vl = insertrightunit (l_LL (ψ, sites[1 ]), 1 ; dual = true )
8052
81- # right side
82- E = @plansor M[1 2 ; 3 ] * Ut[2 ] * ψ. C[sites[end ]][3 ; 4 ] *
83- conj (ψ. C[sites[end ]][1 ; 4 ])
53+ # middle
54+ M = foldl (zip (sites, local_mpo); init = Vl) do v, (site, o)
55+ if o isa Number
56+ return scale! (v * TransferMatrix (ψ. AL[site], ψ. AL[site]), o)
57+ else
58+ return v * TransferMatrix (ψ. AL[site], o, ψ. AL[site])
59+ end
8460 end
8561
86- return E / norm (ψ)^ 2
62+ # right side
63+ E = @plansor removeunit (M, 2 )[1 ; 2 ] * ψ. C[sites[end ]][2 ; 3 ] * conj (ψ. C[sites[end ]][1 ; 3 ])
64+ return E / dot (ψ, ψ)
65+ end
66+
67+ function local_expectation_value1 (ψ:: AbstractMPS , site, O)
68+ E = contract_mpo_expval1 (ψ. AC[site], O, ψ. AC[site])
69+ return E / dot (ψ, ψ)
70+ end
71+ function local_expectation_value2 (ψ:: AbstractMPS , site, O)
72+ AC = ψ. AC[site]
73+ AR = ψ. AR[site + 1 ]
74+ E = contract_mpo_expval2 (AC, AR, O, AC, AR)
75+ return E / dot (ψ, ψ)
8776end
8877
8978# MPOHamiltonian
@@ -93,6 +82,40 @@ function contract_mpo_expval(
9382 )
9483 return @plansor GL[1 2 ; 3 ] * AC[3 7 ; 5 ] * GR[5 8 ; 6 ] * O[2 4 ; 7 8 ] * conj (ACbar[1 4 ; 6 ])
9584end
85+ # generic fallback
86+ function contract_mpo_expval (AC, GL, O, GR, ACbar = AC)
87+ return dot (ACbar, MPO_AC_Hamiltonian (GL, O, GR) * AC)
88+ end
89+
90+ function contract_mpo_expval1 (AC:: MPSTensor , O:: AbstractTensorMap , ACbar:: MPSTensor = AC)
91+ numin (O) == numout (O) == 1 || throw (ArgumentError (" O is not a single-site operator"
92+ return @plansor conj (ACbar[2 3 ; 4 ]) * O[3 ; 1 ] * AC[2 1 ; 4 ]
93+ end
94+ function contract_mpo_expval1 (
95+ AC:: GenericMPSTensor{S, 3} , O:: AbstractTensorMap{<:Any, S} ,
96+ ACbar:: GenericMPSTensor{S, 3} = AC
97+ ) where {S}
98+ numin (O) == numout (O) == 1 || throw (ArgumentError (" O is not a single-site operator"
99+ return @plansor conj (ACbar[2 3 4 ; 5 ]) * O[3 ; 1 ] * AC[2 1 4 ; 5 ]
100+ end
101+
102+ function contract_mpo_expval2 (
103+ A1:: MPSTensor , A2:: MPSTensor , O:: AbstractTensorMap ,
104+ A1bar:: MPSTensor = A1, A2bar:: MPSTensor = A2
105+ )
106+ numin (O) == numout (O) == 2 || throw (ArgumentError (" O is not a two-site operator"
107+ return @plansor conj (A1bar[4 5 ; 6 ]) * conj (A2bar[6 7 ; 8 ]) * O[5 7 ; 2 3 ] * A1[4 2 ; 1 ] *
108+ A2[1 3 ; 8 ]
109+ end
110+ function contract_mpo_expval2 (
111+ A1:: GenericMPSTensor{S, 3} , A2:: GenericMPSTensor{S, 3} ,
112+ O:: AbstractTensorMap{<:Any, S} ,
113+ A1bar:: GenericMPSTensor{S, 3} = A1, A2bar:: GenericMPSTensor{S, 3} = A2
114+ ) where {S}
115+ numin (O) == numout (O) == 2 || throw (ArgumentError (" O is not a two-site operator"
116+ return @plansor conj (A1bar[8 3 4 ; 11 ]) * conj (A2bar[11 12 13 ; 14 ]) * τ[9 6 ; 1 2 ] *
117+ τ[3 4 ; 9 10 ] * A1[8 1 2 ; 5 ] * A2[5 7 13 ; 14 ] * O[10 12 ; 6 7 ]
118+ end
96119
97120function expectation_value (
98121 ψ:: FiniteMPS , H:: FiniteMPOHamiltonian ,
@@ -177,3 +200,12 @@ function expectation_value(
177200 n = norm (ψ. AC[end ])^ 2
178201 return sum (ens) / (n * length (ψ))
179202end
203+
204+ # Density matrices
205+ # ----------------
206+ function expectation_value (ρ:: FiniteMPO , args... )
207+ return expectation_value (convert (FiniteMPS, ρ), args... )
208+ end
209+ function expectation_value (ρ:: InfiniteMPO , args... )
210+ return expectation_value (convert (InfiniteMPS, ρ), args... )
211+ end
0 commit comments