diff --git a/Project.toml b/Project.toml
index 7b664694..acaed477 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "ITensorMPS"
 uuid = "0d1a4710-d33b-49a5-8f18-73bdf49b47e2"
-version = "0.4.0"
+version = "0.4.1"
 authors = ["Matthew Fishman <mfishman@flatironinstitute.org>", "Miles Stoudenmire <mstoudenmire@flatironinstitute.org>"]
 
 [workspace]
diff --git a/src/mps.jl b/src/mps.jl
index 7c768e80..7afe718b 100644
--- a/src/mps.jl
+++ b/src/mps.jl
@@ -907,7 +907,7 @@ function correlation_matrix(
                         Li21 *= oᵢ * dag(psi[pL21])'
                     else
                         sᵢ = siteind(psi, pL21)
-                        Li21 *= prime(dag(si[pL21]), !sᵢ)
+                        Li21 *= prime(dag(psi[pL21]), !sᵢ)
                     end
                     Li21 *= psi[pL21]
                 end
diff --git a/test/base/test_mps.jl b/test/base/test_mps.jl
index 78436cf6..de6947f6 100644
--- a/test/base/test_mps.jl
+++ b/test/base/test_mps.jl
@@ -1004,6 +1004,17 @@ end
 
         @test norm(PM[range] - expect(psi, "S+ * S-"; sites = range)) < 1.0e-8
 
+        # Regression test for https://github.com/ITensor/ITensorMPS.jl/pull/221:
+        # non-contiguous sites with `ishermitian=false` exercise the lower-triangle
+        # contraction path through skipped MPS tensors.
+        non_contiguous = [1, 3, 8]
+        C = correlation_matrix(psi, "Sz", "Sx"; ishermitian = false)
+        Cs =
+            correlation_matrix(psi, "Sz", "Sx"; sites = non_contiguous, ishermitian = false)
+        for (ni, i) in enumerate(non_contiguous), (nj, j) in enumerate(non_contiguous)
+            @test Cs[ni, nj] ≈ C[i, j]
+        end
+
         # With start_site, end_site arguments:
         s = siteinds("S=1/2", 8)
         psi = random_mps(ComplexF64, s; linkdims = m)