Expectation values of local MPO tensors

src/algorithms/expval.jl

@@ -1,22 +1,26 @@
 """
     expectation_value(ψ, O, [environments])
     expectation_value(ψ, inds => O)
+    expectation_value(ψ, (mpo, site => O), [environments])
 
 Compute the expectation value of an operator `O` on a state `ψ`. 
 Optionally, it is possible to make the computations more efficient by also passing in
 previously calculated `environments`.
 
 In general, the operator `O` may consist of an arbitrary MPO `O <: AbstractMPO` that acts
-on all sites, or a local operator `O = inds => operator` acting on a subset of sites. 
-In the latter case, `inds` is a tuple of indices that specify the sites on which the operator
-acts, while the operator is either a `AbstractTensorMap` or a `FiniteMPO`.
+on all sites, a local operator `O = inds => operator` acting on a subset of sites, or a local MPO tensor
+acting on a site within a network whose environment is determined by another MPO `mpo`.
+In the second case, `inds` is a tuple of indices that specify the sites on which the operator
+acts, while the operator is either a `AbstractTensorMap` or a `FiniteMPO`. In the latter case,
+the operator is a `AbstractTensorMap` that acts on the physical space of a single site.
 
 # Arguments
 * `ψ::AbstractMPS` : the state on which to compute the expectation value
-* `O::Union{AbstractMPO,Pair}` : the operator to compute the expectation value of. 
-    This can either be an `AbstractMPO`, or a pair of indices and local operator..
+* `O::Union{AbstractMPO,Pair,AbstractTensorMap}` : the operator to compute the expectation value of. 
+    This can either be an `AbstractMPO`, a pair of indices and local operator, or a local MPO tensor
+    represented as a `AbstractTensorMap`.
 * `environments::AbstractMPSEnvironments` : the environments to use for the calculation. If not given, they will be calculated.
-
+    Depending on the type of `O`, these will be the environments of the operator `O` or the MPO `mpo`.
 # Examples
 ```jldoctest
 julia> ψ = FiniteMPS(ones(Float64, (ℂ^2)^4));
@@ -128,13 +132,24 @@ function expectation_value(
         ψ::InfiniteMPS, H::InfiniteMPOHamiltonian,
         envs::AbstractMPSEnvironments = environments(ψ, H)
     )
-    return sum(1:length(ψ)) do i
+    return sum(1:length(ψ)) do site
         return contract_mpo_expval(
-            ψ.AC[i], envs.GLs[i], H[i][:, 1, 1, end], envs.GRs[i][end]
+            ψ.AC[site], envs.GLs[site], H[site][:, 1, 1, end], envs.GRs[site][end]
         )
     end
 end
 
+# MPO tensor
+#-----------
+function expectation_value(
+        ψ::InfiniteMPS, mpo_pair::Tuple{InfiniteMPO, Pair{Int, <:MPOTensor}},
+        envs::AbstractMPSEnvironments = environments(ψ, first(mpo_pair))
+    )
+    O_mpo, (site, O) = mpo_pair
+    return contract_mpo_expval(ψ.AC[site], envs.GLs[site], O, envs.GRs[site]) /
+        expectation_value(ψ, O_mpo, envs)
+end
+
 # DenseMPO
 # --------
 function expectation_value(ψ::FiniteMPS, mpo::FiniteMPO)

test/algorithms/statmech.jl

@@ -99,4 +99,28 @@ using TensorKit: ℙ
         vals_mps = @constinferred(transfer_spectrum(rho_mps; num_vals))
         @test vals_taylor[1:num_vals] ≈ vals_mps[1:num_vals]
     end
+
+    @testset "2D infinite partition functions with boundary MPS" verbose = true begin
+        beta = 0.5 # ferromagnetic phase
+        f_th = -2.0515856253898357
+        m_th = 0.911319377877496
+        e_th = -1.7455645753125533
+
+        alg = VOMPS(; tol = 1.0e-8, verbosity = 1)
+        O_mpo = classical_ising(; β = beta)
+        ψ₀ = InfiniteMPS(ℂ^2, ℂ^10)
+        ψ, envs = leading_boundary(ψ₀, O_mpo, alg)
+
+        λ = expectation_value(ψ, O_mpo, envs)
+        f = -log(λ) / beta
+        @test f ≈ f_th atol = 1.0e-10
+
+        O, M, E = classical_ising_tensors(beta)
+
+        m = expectation_value(ψ, (O_mpo, 1 => M)) # normalised to give density
+        @test abs(m) ≈ m_th atol = 1.0e-8 # account for spin flip
+
+        e = expectation_value(ψ, (O_mpo, 1 => E))
+        @test e ≈ e_th atol = 1.0e-2
+    end
 end

test/utilities/testsetup.jl

@@ -21,7 +21,7 @@ export force_planar
 export symm_mul_mpo
 export transverse_field_ising, heisenberg_XXX, bilinear_biquadratic_model, XY_model,
     kitaev_model
-export classical_ising, sixvertex
+export classical_ising_tensors, classical_ising, sixvertex
 
 # using TensorOperations
 
@@ -180,24 +180,52 @@ function kitaev_model(
 end
 
 function ising_bond_tensor(β)
-    t = [exp(β) exp(-β); exp(-β) exp(β)]
+    J = 1.0
+    K = β * J
+
+    # Boltzmann weights
+    t = ComplexF64[exp(K) exp(-K); exp(-K) exp(K)]
     r = eigen(t)
     nt = r.vectors * sqrt(Diagonal(r.values)) * r.vectors
     return nt
 end
 
-function classical_ising(; β = log(1 + sqrt(2)) / 2, L = Inf)
-    nt = ising_bond_tensor(β)
+function classical_ising_tensors(beta)
+    nt = ising_bond_tensor(beta)
+    J = 1.0
 
+    # local partition function tensor
     δbulk = zeros(ComplexF64, (2, 2, 2, 2))
     δbulk[1, 1, 1, 1] = 1
     δbulk[2, 2, 2, 2] = 1
     @tensor obulk[-1 -2; -3 -4] := δbulk[1 2; 3 4] * nt[-1; 1] * nt[-2; 2] * nt[-3; 3] *
         nt[-4; 4]
-    Obulk = TensorMap(obulk, ℂ^2 * ℂ^2, ℂ^2 * ℂ^2)
+
+    # magnetization tensor
+    M = copy(δbulk)
+    M[2, 2, 2, 2] *= -1
+    @tensor m[-1 -2; -3 -4] := M[1 2; 3 4] * nt[-1; 1] * nt[-2; 2] * nt[-3; 3] * nt[-4; 4]
+
+    # bond interaction tensor and energy-per-site tensor
+    e = ComplexF64[-J J; J -J] .* nt
+    @tensor e_hor[-1 -2; -3 -4] :=
+        δbulk[1 2; 3 4] * nt[-1; 1] * nt[-2; 2] * nt[-3; 3] * e[-4; 4]
+    @tensor e_vert[-1 -2; -3 -4] :=
+        δbulk[1 2; 3 4] * nt[-1; 1] * nt[-2; 2] * e[-3; 3] * nt[-4; 4]
+    e = e_hor + e_vert
+
+    # fixed tensor map space for all three
+    TMS = ℂ^2 ⊗ ℂ^2 ← ℂ^2 ⊗ ℂ^2
+
+    return TensorMap(obulk, TMS), TensorMap(m, TMS), TensorMap(e, TMS)
+end
+
+function classical_ising(; β = log(1 + sqrt(2)) / 2, L = Inf)
+    Obulk, _, _ = classical_ising_tensors(β)
 
     L == Inf && return InfiniteMPO([Obulk])
 
+    nt = ising_bond_tensor(β)
     δleft = zeros(ComplexF64, (1, 2, 2, 2))
     δleft[1, 1, 1, 1] = 1
     δleft[1, 2, 2, 2] = 1