diff --git a/src/algorithms/expval.jl b/src/algorithms/expval.jl
index e8ac68738..baf32d13f 100644
--- a/src/algorithms/expval.jl
+++ b/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));
@@ -106,12 +110,31 @@ function expectation_value(
         envs::AbstractMPSEnvironments = environments(ψ, H)
     )
     return sum(1:length(ψ)) do i
-        return contract_mpo_expval(
-            ψ.AC[i], envs.GLs[i], H[i][:, 1, 1, end], envs.GRs[i][end]
-        )
+        return _local_expectation_value(ψ, H, envs, i)
     end
 end
 
+function _local_expectation_value(
+        ψ::InfiniteMPS, H::InfiniteMPOHamiltonian,
+        envs::AbstractMPSEnvironments = environments(ψ, H), site::Int = 1
+    )
+    return contract_mpo_expval(
+        ψ.AC[site], envs.GLs[site], H[site][:, 1, 1, end], envs.GRs[site][end]
+    )
+end
+
+# MPO tensor
+#-----------
+function expectation_value(
+        ψ::InfiniteMPS, mpo_pair::Tuple{InfiniteMPO, Pair},
+        envs::AbstractMPSEnvironments = environments(ψ, first(mpo_pair))
+    )
+    site, O = mpo_pair[2]
+    return contract_mpo_expval(
+        ψ.AC[site], envs.GLs[site], O, envs.GRs[site]
+    )
+end
+
 # DenseMPO
 # --------
 function expectation_value(ψ::FiniteMPS, mpo::FiniteMPO)
@@ -128,11 +151,18 @@ function expectation_value(
         envs::MultilineEnvironments = environments(ψ, O)
     )
     return prod(product(1:size(ψ, 1), 1:size(ψ, 2))) do (i, j)
-        GL = envs[i].GLs[j]
-        GR = envs[i].GRs[j]
-        return contract_mpo_expval(ψ.AC[i, j], GL, O[i, j], GR, ψ.AC[i + 1, j])
+        return _local_expectation_value(ψ, O, envs, (i, j))
     end
 end
+function _local_expectation_value(
+        ψ::MultilineMPS, O::MultilineMPO{<:InfiniteMPO},
+        envs::MultilineEnvironments = environments(ψ, O),
+        (i, j)::Tuple{Int, Int} = size(ψ)
+    )
+    GL = envs[i].GLs[j]
+    GR = envs[i].GRs[j]
+    return contract_mpo_expval(ψ.AC[i, j], GL, O[i, j], GR, ψ.AC[i + 1, j])
+end
 function expectation_value(ψ::MultilineMPS, mpo::MultilineMPO, envs...)
     # TODO: fix environments
     return prod(x -> expectation_value(x...), zip(parent(ψ), parent(mpo)))
diff --git a/test/algorithms.jl b/test/algorithms.jl
index b009c121b..64b75a47a 100644
--- a/test/algorithms.jl
+++ b/test/algorithms.jl
@@ -462,6 +462,34 @@ module TestAlgorithms
         end
     end
 
+    @testset "2D infinite partition functions" 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 = verbosity_conv)
+        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)
+
+        denom = expectation_value(ψ, O_mpo, envs)
+        denom2 = expectation_value(ψ, (O_mpo, 1 => O), envs)
+        @test denom ≈ denom2 atol = 1.0e-10
+
+        m = expectation_value(ψ, (O_mpo, 1 => M)) / denom
+        @test abs(m) ≈ m_th atol = 1.0e-8 # account for spin flip
+
+        e = expectation_value(ψ, (O_mpo, 1 => E)) / denom
+        @test e ≈ e_th atol = 1.0e-2
+    end
+
     @testset "excitations" verbose = true begin
         @testset "infinite (ham)" begin
             H = repeat(force_planar(heisenberg_XXX()), 2)
diff --git a/test/setup.jl b/test/setup.jl
index a6b07bbc0..8c5cf92e2 100644
--- a/test/setup.jl
+++ b/test/setup.jl
@@ -18,7 +18,7 @@ export force_planar
 export symm_mul_mpo
 export transverse_field_ising, heisenberg_XXX, bilinear_biquadratic_model, XY_model,
     kitaev_model
-export classical_ising, finite_classical_ising, sixvertex
+export classical_ising_tensors, classical_ising, sixvertex
 
 # using TensorOperations
 
@@ -294,24 +294,52 @@ function kitaev_model(; t = 1.0, mu = 1.0, Delta = 1.0, L = Inf)
 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