From b6d51d28ba29513a49425db5cddfb63961e3d1dd Mon Sep 17 00:00:00 2001
From: Ben Corbett <corbett5@llnl.gov>
Date: Sun, 26 Apr 2026 22:13:26 -0600
Subject: [PATCH 01/16] Made a plan

---
 plan.md | 518 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 518 insertions(+)
 create mode 100644 plan.md

diff --git a/plan.md b/plan.md
new file mode 100644
index 0000000..b4ce6f3
--- /dev/null
+++ b/plan.md
@@ -0,0 +1,518 @@
+# Symbolic Vertex-Cover MPO Construction
+
+## Summary
+
+Add a reusable symbolic MPO construction path for integer-labeled `OpIDSum`s.
+`MPO_symbolic` should run the vertex-cover MPO construction once, store the
+resulting tensor entries as integer linear combinations of coefficient
+variables, and allow repeated numeric instantiation with different coefficient
+vectors.
+
+The primary API should follow the package's current argument order:
+
+```julia
+sym = MPO_symbolic(os, sites; basis_op_cache_vec=nothing,
+                   check_for_errors=true,
+                   combine_qn_sectors=false,
+                   output_level=0)
+
+H = instantiate_MPO(sym, coefficients; splitblocks=true, checkflux=false)
+instantiate_MPO!(H, sym, coefficients; checkflux=false)
+H2 = instantiate_MPO(H, sym, coefficients; checkflux=false)
+```
+
+## Public API and V1 Semantics
+
+- Export `SymbolicMPO` and `MPO_symbolic`.
+- Support only `OpIDSum{N,<:Integer,Ti}` input in V1.
+- Symbolic construction is vertex-cover construction. It should not expose an
+  `alg` keyword, and there is no QR symbolic path to support.
+- Do not add symbolic `OpSum` conversion in V1.
+- Preserve numeric `MPO_new` behavior exactly.
+- User-provided integer labels are variable labels and must always be greater
+  than zero.
+- Convert each user label `k` to internal symbolic id `k + 1`, reserving
+  internal id `1` for the constant one introduced by the vertex-cover algorithm.
+- Carry signs by negating the internal symbolic id. Negative ids can arise from
+  fermionic sorting signs or sign-only basis rewrites; they are not user-facing
+  labels.
+- During numeric substitution:
+  - internal id `+1` evaluates to `+1.0`,
+  - internal id `-1` evaluates to `-1.0`,
+  - internal id `k > 1` evaluates to `coefficients[k - 1]`,
+  - internal id `k < -1` evaluates to `-coefficients[-k - 1]`.
+- The coefficient vector is indexed by the original positive user labels. The
+  reserved internal constant id is transparent to the user, so there is no
+  reserved slot in `coefficients`.
+- Reject `OpSum` input, unsupported basis rewrite factors, and missing
+  coefficients with clear `ArgumentError`s.
+
+## Internal Representation
+
+### Symbolic local matrix terms
+
+Do not model symbolic MPO entries as scalar coefficients. A symbolic block entry
+is a symbolic weighted sum of local operator matrices:
+
+```julia
+const SymbolicLocalMatrix{Ti<:Integer} = Vector{Tuple{Int,Ti}}
+```
+
+Each tuple is `(signed_weight_id, signed_local_op_id)`.
+
+- `signed_weight_id` is the signed internal coefficient id:
+  - `+1` and `-1` are the constants `+1` and `-1`,
+  - `k > 1` maps to user assignment `coefficients[k - 1]`,
+  - `k < -1` maps to `-coefficients[-k - 1]`.
+- `abs(signed_local_op_id)::Ti` uses the same integer type as
+  `OpIDSum{N,C,Ti}` and identifies the local operator being scaled by that
+  signed coefficient id.
+- `signed_local_op_id < 0` means the local operator should be emitted with the
+  Jordan-Wigner string applied. `signed_local_op_id > 0` means no
+  Jordan-Wigner string is applied.
+- A term represents
+  `substitute(signed_weight_id, coefficients) * local_matrix`, where
+  `local_matrix` is `op_cache_vec[n][abs(signed_local_op_id)].matrix`, with the
+  Jordan-Wigner sign pattern applied when `signed_local_op_id < 0`.
+
+Required helpers:
+
+- `_weight_value(signed_weight_id, coefficients)` substitutes one signed
+  internal id using the rules above.
+- `_normalize_symbolic_local_matrix!(terms)` canonicalizes ordering and removes
+  canceling pairs such as `(k, op)` with `(-k, op)`. It must not collapse
+  repeated identical pairs into one entry, since multiplicity represents
+  repeated contributions such as `2 * coefficients[k - 1] * op`.
+- `_scale_symbolic_weight(signed_weight_id, factor)` accepts only exact factors
+  `+1` and `-1` in V1, flipping the signed id for `-1`.
+- `_evaluate_symbolic_local_matrix(terms, coefficients, op_cache)` builds the
+  dense numeric local matrix by summing cached local operator matrices, applying
+  Jordan-Wigner strings for negative local operator ids, and substituting
+  weights.
+
+### Symbolic Block Storage
+
+Do not reuse `BlockSparseMatrix{C}` for symbolic block entries, since that alias
+stores dense `Matrix{C}` blocks. Symbolic construction should use:
+
+```julia
+const SymbolicBlockSparseMatrix{Ti<:Integer} =
+  Vector{Dictionary{Int,SymbolicLocalMatrix{Ti}}}
+```
+
+This mirrors the existing `right_link => left_link => block` layout, but each
+block is a symbolic local matrix rather than an already-expanded dense matrix.
+Numeric instantiation converts these symbolic blocks into ordinary
+`BlockSparseMatrix{T}` values immediately before constructing `ITensor`s.
+
+## Implementation Checklist
+
+### 1. Add and test symbolic local matrix operations
+
+- Add `SymbolicLocalMatrix{Ti}` and `SymbolicBlockSparseMatrix{Ti}` aliases or
+  small wrapper structs.
+- Implement helpers for appending terms, normalizing terms, applying `+1` or
+  `-1` rewrite factors, computing `max_user_label`, and substituting signed
+  internal ids.
+- Implement `_evaluate_symbolic_local_matrix` to produce a dense numeric matrix
+  for one site from `(signed_weight_id, signed_local_op_id)` terms and the site's
+  operator cache.
+- Add focused tests for:
+  - duplicate `(weight_id, signed_local_op_id)` terms preserving multiplicity,
+  - zero normalization,
+  - `+1` and `-1` constant substitution,
+  - positive user label to internal id mapping,
+  - signed internal id evaluation,
+  - missing coefficient rejection,
+  - non-`±1` rewrite-factor rejection.
+
+Tests for this step:
+
+- Unit-test `_weight_value` directly:
+  - `+1` returns `1`,
+  - `-1` returns `-1`,
+  - `+3` returns `coefficients[2]`,
+  - `-3` returns `-coefficients[2]`,
+  - `+4` throws when `length(coefficients) == 2`.
+- Unit-test symbolic local matrix normalization:
+  - `[(3, 2), (-3, 2)]` normalizes to an empty term list,
+  - `[(3, 2), (3, 2)]` keeps two entries and evaluates to twice the same local
+    matrix contribution,
+  - terms are sorted deterministically so equality checks are stable.
+- Unit-test Jordan-Wigner encoding:
+  - `(3, -2)` evaluates to `coefficients[2]` times operator `2` with the
+    Jordan-Wigner sign pattern,
+  - `(3, 2)` evaluates to the same operator without the Jordan-Wigner sign
+    pattern.
+
+### 2. Convert the integer-labeled `OpIDSum` into internal-id form
+
+- Add an internal conversion helper:
+
+```julia
+internalize_symbolic_ids!(os::OpIDSum{N,C,Ti}) where {N,C<:Integer,Ti}
+```
+
+- This helper may mutate the provided `OpIDSum`, matching the existing MPO
+  construction behavior where preprocessing, sorting, and duplicate compaction
+  are destructive.
+- Map each stored scalar label to a signed internal id in place.
+- Reject scalar label `0`.
+- Map stored scalar label `k > 0` to internal id `k + 1`.
+- Map stored scalar label `k < 0` to internal id `-(abs(k) + 1)`, treating it
+  as the sign-carrying form of user label `abs(k)`.
+- Keep the same `op_cache_vec`, `abs_tol`, and `modify!`.
+- Add tests proving the in-place conversion maps labels correctly and rejects
+  zero labels.
+
+Tests for this step:
+
+- Build a small `OpIDSum{2,Int,Int}` with scalar labels `1`, `2`, and `-3`;
+  after `internalize_symbolic_ids!`, assert that the stored labels are `2`, `3`,
+  and `-4`.
+- Add a term with label `0` and assert `internalize_symbolic_ids!` throws an
+  `ArgumentError`.
+- Assert the helper is intentionally mutating by checking that the same `os`
+  object has updated scalar storage after conversion.
+- Assert `op_cache_vec`, `abs_tol`, and `modify!` are still the original objects
+  or equivalent values after conversion.
+
+### 3. Make basis rewrite compatible with internal signed ids
+
+- Reuse `prepare_opID_sum!(os, basis_op_cache_vec)` after the in-place
+  internal-id conversion.
+- Ensure `rewrite_in_operator_basis!` can multiply internal signed ids by basis
+  rewrite factors only when the factor is exactly `+1` or `-1`.
+- Reject `im`, `2im`, `0.5`, `2`, and other unsupported factors with a clear
+  `ArgumentError`.
+- Add a sign-only rewrite regression where a same-site product produces `-1`
+  and verifies the internal id sign flips while the user label is preserved.
+
+Tests for this step:
+
+- Construct a small operator-cache fixture where a same-site product rewrites
+  into a basis operator with factor `-1`; after preparation, verify internal id
+  `k + 1` becomes `-(k + 1)` and the basis operator id is updated.
+- Construct a fixture with rewrite factor `+1`; verify the internal id sign is
+  unchanged.
+- Reuse the existing `X * Y -> im * Z` style case or an equivalent fixture to
+  assert complex factors throw.
+- Add separate fixtures for `0.5` and `2` factors and assert both throw, since
+  V1 accepts only exact `±1`.
+
+### 4. Build a symbolic graph without adding identifier integers
+
+- Numeric `MPOGraph(os::OpIDSum)` must keep its current duplicate-term
+  compaction and tolerance behavior.
+- Add a symbolic graph-construction path for the internal-id `OpIDSum`.
+- Duplicate operator tuples must combine correctly even when their scalar labels
+  are different symbolic variables. Combining duplicates should preserve all
+  signed ids, for example by storing a vector of signed internal ids as the edge
+  weight or by expanding the duplicates into distinct edge entries before local
+  block construction.
+- The graph edge weight type should remain a signed internal id or a collection
+  of signed ids, not a symbolic local matrix term.
+- Add a test where duplicate operator tuples with different user labels both
+  contribute to the same symbolic block and instantiate correctly.
+
+Tests for this step:
+
+- Build an `OpIDSum{2,Int,Int}` with two identical operator tuples labeled `1`
+  and `2`; after symbolic graph/block construction, assert the resulting
+  symbolic local matrix contains both signed ids, not their numeric sum.
+- Add two identical operator tuples with the same label `1`; instantiate with
+  `coefficients[1] = a` and verify the numeric MPO contribution is `2a`.
+- Add two identical operator tuples whose labels become opposite signed internal
+  ids after preprocessing; verify the symbolic local matrix cancels or
+  instantiates to zero for that contribution.
+- Run the same duplicate-label test through `MPO_symbolic` and compare the
+  instantiated MPO against numeric `MPO_new` built from the corresponding
+  numeric coefficients.
+
+### 5. Build `SymbolicMPO`
+
+Add:
+
+```julia
+struct SymbolicMPO{Ti<:Integer}
+  offsets::Vector{Vector{Int}}
+  block_sparse_matrices::Vector{Vector{SymbolicBlockSparseMatrix{Ti}}}
+  sites::Vector{<:Index}
+  llinks::Vector{Index}
+  op_cache_vec::OpCacheVec
+  max_user_label::Int
+end
+```
+
+If Julia requires concrete field types for `sites`, parameterize the site vector
+type instead of using an abstract field.
+
+`MPO_symbolic(os, sites; kwargs...)` should:
+
+- validate `os` has an integer coefficient type,
+- reject any supplied `alg` keyword with an `ArgumentError` explaining that
+  symbolic construction is always vertex-cover based,
+- convert `basis_op_cache_vec` with `to_OpCacheVec(sites, basis_op_cache_vec)`,
+- convert `os` to internal-id form in place,
+- run `prepare_opID_sum!` and optionally `check_os_for_errors`,
+- build the symbolic graph without numerically adding identifier ids,
+- initialize the dummy left link exactly as `MPO_new` does,
+- run a VC-only symbolic construction driver that mirrors the existing
+  vertex-cover path but writes `SymbolicBlockSparseMatrix{Ti}` blocks,
+- when adding a local contribution, append
+  `(signed_weight_id, signed_local_op_id)` to the symbolic block instead of
+  expanding the dense local operator matrix,
+- store offsets, symbolic block matrices, sites, links, `op_cache_vec`, and the
+  maximum positive user label needed for substitution.
+
+Tests for this step:
+
+- Construct a small transverse-field Ising `OpIDSum{2,Int,Int}` and assert
+  `MPO_symbolic(os, sites)` returns a `SymbolicMPO` with:
+  - `length(offsets) == length(sites)`,
+  - `length(block_sparse_matrices) == length(sites)`,
+  - `length(llinks) == length(sites) + 1`,
+  - `max_user_label` equal to the largest user label in the input.
+- Assert `MPO_symbolic(os, sites; alg="QR")` throws an `ArgumentError`.
+- Assert `MPO_symbolic(::OpSum, sites)` throws an `ArgumentError` or has no
+  method, whichever API choice is implemented.
+- For a fermionic fixture, inspect at least one symbolic local matrix term and
+  assert negative `signed_local_op_id` is used where `needs_JW_string` is true.
+
+### 6. Instantiate a fresh numeric MPO
+
+Add:
+
+```julia
+instantiate_MPO(sym::SymbolicMPO, coefficients; splitblocks=true, checkflux=false)
+```
+
+It should:
+
+- validate the coefficient vector length against `sym.max_user_label`,
+- evaluate symbolic block entries into numeric block matrices using
+  `sym.op_cache_vec`,
+- preserve the symbolic structural pattern even when evaluated values are zero,
+- call the existing tensor assembly path with the requested `splitblocks` value,
+- contract away the dummy boundary links the same way as the low-level
+  `instantiate_MPO(offsets, block_sparse_matrices, sites, llinks; ...)`,
+- optionally run `ITensors.checkflux`.
+
+For QN tensors, structural block discovery must be based on symbolic local
+matrix terms before numeric substitution, so a coefficient vector containing
+zeros does not change the block structure or link layout.
+
+Tests for this step:
+
+- Instantiate a small symbolic MPO with two different coefficient vectors and
+  compare each result against numeric `MPO_new` built from an equivalent numeric
+  `OpIDSum`.
+- Instantiate with `splitblocks=false` and `splitblocks=true` for a QN-conserving
+  fixture and compare both against numeric `MPO_new`.
+- Instantiate with a coefficient vector containing zeros and assert the link
+  dimensions match a nonzero-coefficient instantiation from the same
+  `SymbolicMPO`.
+- Pass a too-short coefficient vector and assert an `ArgumentError`.
+- Run with `checkflux=true` on a QN fixture and assert it succeeds.
+
+### 7. Add in-place and template-assisted instantiation
+
+Add:
+
+```julia
+instantiate_MPO!(H::MPO, sym::SymbolicMPO, coefficients; checkflux=false)
+instantiate_MPO(H_template::MPO, sym::SymbolicMPO, coefficients; checkflux=false)
+```
+
+These overloads must not accept a `splitblocks` keyword. A preconstructed MPO
+already fixes the block layout, so evaluation must use the layout of `H` or
+`H_template`.
+
+`instantiate_MPO!` should:
+
+- check `length(H) == length(sym.sites)`,
+- check site indices are compatible with `sym.sites`,
+- check link dimensions match `sym.llinks`,
+- overwrite each tensor in `H` with the evaluated tensor for that site,
+- preserve the existing `MPO` object identity.
+
+`instantiate_MPO(H_template, sym, coefficients)` should:
+
+- reject incompatible templates with `ArgumentError`,
+- copy or reuse a compatible template,
+- call `instantiate_MPO!`,
+- return the instantiated MPO.
+
+It is acceptable for V1 to overwrite tensors rather than doing lower-level
+storage mutation, as long as the public behavior and link structure are stable.
+
+Tests for this step:
+
+- Build `H = instantiate_MPO(sym, coeffs1; splitblocks=true)`, then call
+  `instantiate_MPO!(H, sym, coeffs2)` and compare `H` against
+  `instantiate_MPO(sym, coeffs2; splitblocks=true)`.
+- Assert `instantiate_MPO!(H, sym, coeffs; splitblocks=false)` throws a
+  `MethodError` or `ArgumentError`, since this overload must not accept
+  `splitblocks`.
+- Call `instantiate_MPO(H, sym, coeffs2)` with a compatible template and compare
+  against fresh instantiation using the same layout.
+- Build an incompatible template with different sites or link dimensions and
+  assert `instantiate_MPO(H_bad, sym, coeffs)` throws an `ArgumentError`.
+
+### 8. Add symbolic MPO tests
+
+Create `test/symbolic-mpo-tests.jl` and include it from `test/runtests.jl`.
+
+The file should use named `@testset`s matching the implementation steps, so a
+failure points directly to the layer that broke:
+
+- `"symbolic local matrix terms"`:
+  - normalization,
+  - duplicate `(weight_id, signed_local_op_id)` multiplicity,
+  - constant weight ids,
+  - negative local operator id applies the Jordan-Wigner string,
+  - positive user label to internal id mapping,
+  - signed internal ids,
+  - missing coefficient error,
+  - unsupported rewrite-factor error.
+- `"internal symbolic ids"`:
+  - in-place user-label to internal-id conversion,
+  - zero label rejection,
+  - sign-preserving conversion for labels made negative by preprocessing.
+- `"basis rewrite"`:
+  - `+1` rewrite factor preserves id sign,
+  - `-1` rewrite factor flips id sign,
+  - complex, fractional, and integer factors other than `±1` throw.
+- `"symbolic graph duplicate terms"`:
+  - different labels on the same operator tuple are both preserved,
+  - repeated identical labels preserve multiplicity,
+  - opposite signed ids cancel or instantiate to zero.
+- `"fresh instantiation"`:
+  - build `OpIDSum{2,Int,Int}`,
+  - construct `sym = MPO_symbolic(os, sites)`,
+  - instantiate with several random coefficient vectors indexed directly by the
+    positive user labels, passing `splitblocks` to fresh `instantiate_MPO`,
+  - compare against numeric `MPO_new` built with the same coefficients.
+- `"public API rejections"`:
+  - any supplied `alg` keyword,
+  - `OpSum` input,
+  - missing coefficient,
+  - incompatible preconstructed MPO templates.
+- `"in-place instantiation"`:
+  - instantiate once,
+  - update with a second coefficient vector through `instantiate_MPO!`,
+  - verify `instantiate_MPO!` does not accept a `splitblocks` keyword,
+  - compare to fresh instantiation,
+  - check link dimensions remain unchanged.
+- `"template-assisted instantiation"`:
+  - instantiate from a compatible template,
+  - compare to fresh instantiation,
+  - reject an incompatible template.
+
+Tests for this step:
+
+- `julia --project test/symbolic-mpo-tests.jl` runs all new symbolic testsets
+  by itself.
+- `julia --project test/runtests.jl` includes `test/symbolic-mpo-tests.jl` once
+  and still runs the existing test files.
+- Existing `test/ops-tests.jl` still passes, proving the basis-rewrite changes
+  did not regress the numeric `OpIDSum` path.
+
+### 9. Update electronic-structure example
+
+Add a compact symbolic section to `examples/electronic-structure.jl`.
+
+Implementation outline:
+
+- Build a small integer-labeled `OpIDSum` in the same term order as the numeric
+  electronic-structure `OpIDSum`.
+- Keep a `coefficients` vector indexed directly by positive user labels.
+- Whenever user label `k` is added to the symbolic `OpIDSum`, store the actual
+  numeric coefficient at `coefficients[k]`.
+- Construct `sym = MPO_symbolic(
+  os_symbolic, sites; basis_op_cache_vec=os_symbolic.op_cache_vec, check_for_errors=false
+  )`.
+- Instantiate twice with two coefficient vectors, passing `splitblocks=true` to
+  fresh `instantiate_MPO`.
+- For a small `N`, compare one symbolic instantiation to the existing numeric
+  `MPO_new` result.
+- Keep the example section short so the existing benchmark narrative remains
+  readable.
+
+Tests for this step:
+
+- Run the example with a small `N` setting first and assert the symbolic
+  instantiation matches the numeric `MPO_new` result using the existing MPO
+  comparison helper or an equivalent relative-difference check.
+- Run the final example command:
+
+```bash
+julia --project -t8 examples/electronic-structure.jl
+```
+
+- Confirm the example still reports the existing numeric construction metrics
+  and additionally prints or checks the symbolic construction section.
+- Confirm changing the second coefficient vector changes the instantiated MPO
+  without rebuilding the symbolic structure.
+
+### 10. Update docs
+
+- Add docstrings for `SymbolicMPO`, `MPO_symbolic`, and the new
+  `instantiate_MPO` / `instantiate_MPO!` overloads.
+- Update `docs/src/documentation/MPO_new.md` to include the symbolic API in the
+  `@docs` block.
+- Add a short prose section explaining:
+  - integer labels in `OpIDSum`,
+  - positive user labels and transparent internal remapping,
+  - repeated instantiation,
+  - why symbolic construction is vertex-cover based and has no QR option,
+  - `±1` rewrite-factor limitation,
+  - `OpIDSum`-only limitation.
+
+Tests for this step:
+
+- Verify `julia --project -e 'using ITensorMPOConstruction'` loads cleanly after
+  adding exports and docstrings.
+- Verify the docs build:
+
+```bash
+julia --project=docs --color=yes docs/make.jl
+```
+
+- Check the generated docs page includes `SymbolicMPO`, `MPO_symbolic`,
+  `instantiate_MPO(sym, ...)`, and `instantiate_MPO!(H, sym, ...)`.
+
+## Verification Checklist
+
+Run focused checks first:
+
+```bash
+julia --project -e 'using ITensorMPOConstruction'
+julia --project test/symbolic-mpo-tests.jl
+julia --project test/ops-tests.jl
+julia --project test/test-MPOConstruction.jl
+```
+
+Then run broader checks:
+
+```bash
+julia --project test/runtests.jl
+julia --project=docs --color=yes docs/make.jl
+```
+
+Because `examples/electronic-structure.jl` is modified, also run:
+
+```bash
+julia --project -t8 examples/electronic-structure.jl
+```
+
+## Acceptance Criteria
+
+- Existing numeric `MPO_new` tests pass unchanged.
+- `MPO_symbolic` may mutate its input `OpIDSum`, consistently with existing MPO
+  construction preprocessing.
+- Repeated symbolic instantiation matches fresh numeric MPO construction.
+- In-place symbolic instantiation updates an existing MPO without changing link
+  dimensions.
+- Symbolic construction rejects unsupported input cases with clear errors.
+- Docs build successfully and include the new public API.

From f2dff1ddc9c3ce4c753a67aba6f96842acbaf6a4 Mon Sep 17 00:00:00 2001
From: Ben Corbett <corbett5@llnl.gov>
Date: Sun, 26 Apr 2026 22:33:53 -0600
Subject: [PATCH 02/16] Step 1

---
 src/ITensorMPOConstruction.jl       |   1 +
 src/OpIDSum.jl                      |   2 -
 src/SymbolicLocalMatrix.jl          | 185 ++++++++++++++++++++++++++++
 test/runtests.jl                    |   1 +
 test/symbolic-local-matrix-tests.jl |  79 ++++++++++++
 5 files changed, 266 insertions(+), 2 deletions(-)
 create mode 100644 src/SymbolicLocalMatrix.jl
 create mode 100644 test/symbolic-local-matrix-tests.jl

diff --git a/src/ITensorMPOConstruction.jl b/src/ITensorMPOConstruction.jl
index c81a4f9..b3b37bb 100644
--- a/src/ITensorMPOConstruction.jl
+++ b/src/ITensorMPOConstruction.jl
@@ -18,6 +18,7 @@ using ThreadsX
 #
 include("time_if.jl")
 include("OpIDSum.jl")
+include("SymbolicLocalMatrix.jl")
 include("ops.jl")
 include("BipartiteGraph.jl")
 include("connected-components.jl")
diff --git a/src/OpIDSum.jl b/src/OpIDSum.jl
index e02eef5..fc4c36e 100644
--- a/src/OpIDSum.jl
+++ b/src/OpIDSum.jl
@@ -39,8 +39,6 @@ function OpInfo(local_op::Op, site::Index)
   )
 end
 
-# TODO: Create a symbolic sum.
-
 """
     OpCacheVec
 
diff --git a/src/SymbolicLocalMatrix.jl b/src/SymbolicLocalMatrix.jl
new file mode 100644
index 0000000..aecf146
--- /dev/null
+++ b/src/SymbolicLocalMatrix.jl
@@ -0,0 +1,185 @@
+"""
+    SymbolicLocalMatrix{Ti}
+
+Symbolic weighted sum of cached local operators for one MPO block entry.
+
+Each term is `(signed_weight_id, signed_local_op_id)`. The signed weight id is
+an internal symbolic coefficient id, where `+1` and `-1` represent constants,
+and `+k` or `-k` for `k > 1` map to user coefficient `coefficients[k - 1]` with
+the stored sign. The absolute value of `signed_local_op_id` is the local
+operator id; negative local operator ids indicate that a Jordan-Wigner string
+should be applied while evaluating the cached operator matrix.
+"""
+const SymbolicLocalMatrix{Ti<:Integer} = Vector{Tuple{Int,Ti}}
+
+"""
+    SymbolicBlockSparseMatrix{Ti}
+
+Block-sparse symbolic storage mirroring `BlockSparseMatrix`, but with each
+block entry represented as a [`SymbolicLocalMatrix`](@ref).
+"""
+const SymbolicBlockSparseMatrix{Ti<:Integer} =
+  Vector{Dictionary{Int,SymbolicLocalMatrix{Ti}}}
+
+function _check_symbolic_weight_id(signed_weight_id::Integer)::Nothing
+  signed_weight_id == 0 && throw(
+    ArgumentError(
+      "Symbolic weight id 0 is invalid. Internal id 1 is reserved for the constant one.",
+    ),
+  )
+  return nothing
+end
+
+function _check_symbolic_local_op_id(signed_local_op_id::Integer)::Nothing
+  signed_local_op_id == 0 && throw(
+    ArgumentError("Symbolic local operator id 0 is invalid. Operator ids are one-based."),
+  )
+  return nothing
+end
+
+function _internal_symbolic_id(user_label::Integer)::Int
+  user_label == 0 && throw(ArgumentError("Symbolic coefficient labels must be nonzero."))
+  return user_label > 0 ? user_label + 1 : -(abs(user_label) + 1)
+end
+
+function _max_user_label(signed_weight_id::Integer)::Int
+  _check_symbolic_weight_id(signed_weight_id)
+  return abs(signed_weight_id) - 1
+end
+
+function _max_user_label(terms::SymbolicLocalMatrix)::Int
+  max_user_label = 0
+  for (signed_weight_id, _) in terms
+    max_user_label = max(max_user_label, _max_user_label(signed_weight_id))
+  end
+  return max_user_label
+end
+
+function _append_symbolic_local_matrix_term!(
+  terms::SymbolicLocalMatrix{Ti},
+  signed_weight_id::Integer,
+  signed_local_op_id::Integer,
+)::SymbolicLocalMatrix{Ti} where {Ti}
+  _check_symbolic_weight_id(signed_weight_id)
+  _check_symbolic_local_op_id(signed_local_op_id)
+  push!(terms, (Int(signed_weight_id), Ti(signed_local_op_id)))
+  return terms
+end
+
+function _normalize_symbolic_local_matrix!(
+  terms::SymbolicLocalMatrix{Ti}
+)::SymbolicLocalMatrix{Ti} where {Ti}
+  for (signed_weight_id, signed_local_op_id) in terms
+    _check_symbolic_weight_id(signed_weight_id)
+    _check_symbolic_local_op_id(signed_local_op_id)
+  end
+
+  sort!(terms; by=term -> (term[2], abs(term[1]), term[1]))
+
+  read_index = 1
+  write_index = 1
+  while read_index <= length(terms)
+    signed_local_op_id = terms[read_index][2]
+    abs_weight_id = abs(terms[read_index][1])
+    num_positive = 0
+    num_negative = 0
+
+    while read_index <= length(terms) &&
+        terms[read_index][2] == signed_local_op_id &&
+        abs(terms[read_index][1]) == abs_weight_id
+      if terms[read_index][1] > 0
+        num_positive += 1
+      else
+        num_negative += 1
+      end
+      read_index += 1
+    end
+
+    multiplicity = num_positive - num_negative
+    signed_weight_id = multiplicity > 0 ? abs_weight_id : -abs_weight_id
+    for _ in 1:abs(multiplicity)
+      terms[write_index] = (signed_weight_id, signed_local_op_id)
+      write_index += 1
+    end
+  end
+
+  resize!(terms, write_index - 1)
+  return terms
+end
+
+function _scale_symbolic_weight(signed_weight_id::Integer, factor)::Int
+  _check_symbolic_weight_id(signed_weight_id)
+  isone(factor) && return Int(signed_weight_id)
+  factor == -one(factor) && return -Int(signed_weight_id)
+  throw(
+    ArgumentError(
+      "Unsupported symbolic rewrite factor $factor. " *
+      "Symbolic MPO construction supports only exact +1 and -1 factors.",
+    ),
+  )
+end
+
+function _weight_value(signed_weight_id::Integer, coefficients::AbstractVector)
+  _check_symbolic_weight_id(signed_weight_id)
+  abs_weight_id = abs(signed_weight_id)
+  abs_weight_id == 1 && return signed_weight_id > 0 ? 1.0 : -1.0
+
+  user_label = abs_weight_id - 1
+  user_label <= length(coefficients) || throw(
+    ArgumentError(
+      "Missing coefficient for symbolic label $user_label. " *
+      "Received $(length(coefficients)) coefficients.",
+    ),
+  )
+
+  value = coefficients[user_label]
+  return signed_weight_id > 0 ? value : -value
+end
+
+function _symbolic_local_matrix_eltype(
+  coefficients::AbstractVector, op_cache::Vector{OpInfo}
+)::Type
+  val_type = promote_type(Float64, eltype(coefficients))
+  for op_info in op_cache
+    val_type = promote_type(val_type, eltype(op_info.matrix))
+  end
+  return val_type
+end
+
+function _evaluate_symbolic_local_matrix(
+  terms::SymbolicLocalMatrix, coefficients::AbstractVector, op_cache::Vector{OpInfo}
+)::Matrix
+  isempty(op_cache) && throw(
+    ArgumentError("Cannot evaluate a symbolic local matrix with an empty operator cache."),
+  )
+
+  matrix_size = size(op_cache[1].matrix)
+  result = zeros(_symbolic_local_matrix_eltype(coefficients, op_cache), matrix_size)
+
+  for (signed_weight_id, signed_local_op_id) in terms
+    _check_symbolic_weight_id(signed_weight_id)
+    _check_symbolic_local_op_id(signed_local_op_id)
+
+    local_op_id = abs(Int(signed_local_op_id))
+    local_op_id <= length(op_cache) || throw(
+      ArgumentError(
+        "Symbolic local operator id $local_op_id is not available in a cache with " *
+        "$(length(op_cache)) operators.",
+      ),
+    )
+
+    local_op = op_cache[local_op_id].matrix
+    size(local_op) == matrix_size || throw(
+      ArgumentError(
+        "Operators in a symbolic local matrix must have matching matrix sizes. " *
+        "Expected $matrix_size, got $(size(local_op)).",
+      ),
+    )
+
+    add_to_local_matrix!(
+      result, _weight_value(signed_weight_id, coefficients), local_op, signed_local_op_id < 0
+    )
+  end
+
+  return result
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 3541f69..fd80071 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,5 +1,6 @@
 include("tensor-construction-tests.jl")
 include("ops-tests.jl")
+include("symbolic-local-matrix-tests.jl")
 include("large-graph-tests.jl")
 include("test-MPOConstruction.jl")
 include("mponew_h1_accuracy.jl")
diff --git a/test/symbolic-local-matrix-tests.jl b/test/symbolic-local-matrix-tests.jl
new file mode 100644
index 0000000..56d76dd
--- /dev/null
+++ b/test/symbolic-local-matrix-tests.jl
@@ -0,0 +1,79 @@
+using ITensorMPOConstruction
+using ITensorMPOConstruction:
+  SymbolicLocalMatrix,
+  _append_symbolic_local_matrix_term!,
+  _evaluate_symbolic_local_matrix,
+  _internal_symbolic_id,
+  _max_user_label,
+  _normalize_symbolic_local_matrix!,
+  _scale_symbolic_weight,
+  _weight_value
+using ITensors
+using ITensorMPS
+using Test
+
+function test_symbolic_local_matrix_operations()::Nothing
+  coefficients = [10.0, 20.0]
+
+  @test _weight_value(1, coefficients) == 1.0
+  @test _weight_value(-1, coefficients) == -1.0
+  @test _weight_value(3, coefficients) == coefficients[2]
+  @test _weight_value(-3, coefficients) == -coefficients[2]
+  @test_throws ArgumentError _weight_value(4, coefficients)
+
+  @test _internal_symbolic_id(1) == 2
+  @test _internal_symbolic_id(3) == 4
+  @test_throws ArgumentError _internal_symbolic_id(0)
+
+  @test _scale_symbolic_weight(3, 1) == 3
+  @test _scale_symbolic_weight(3, -1) == -3
+  @test _scale_symbolic_weight(-3, -1) == 3
+  @test_throws ArgumentError _scale_symbolic_weight(3, 2)
+  @test_throws ArgumentError _scale_symbolic_weight(3, im)
+
+  terms = SymbolicLocalMatrix{Int}()
+  _append_symbolic_local_matrix_term!(terms, 3, -2)
+  @test terms == [(3, -2)]
+  @test _max_user_label(SymbolicLocalMatrix{Int}([(1, 1), (-3, 2), (5, -2)])) == 4
+
+  terms = SymbolicLocalMatrix{Int}([(3, 2), (-3, 2)])
+  _normalize_symbolic_local_matrix!(terms)
+  @test isempty(terms)
+
+  terms = SymbolicLocalMatrix{Int}([(3, 2), (3, 2)])
+  _normalize_symbolic_local_matrix!(terms)
+  @test terms == [(3, 2), (3, 2)]
+
+  terms = SymbolicLocalMatrix{Int}([(5, 3), (4, 2), (3, 2)])
+  _normalize_symbolic_local_matrix!(terms)
+  @test terms == [(3, 2), (4, 2), (5, 3)]
+
+  qubit_sites = siteinds("Qubit", 1)
+  qubit_op_cache_vec = to_OpCacheVec(qubit_sites, [["I", "X"]])
+  duplicate_terms = SymbolicLocalMatrix{Int}([(3, 2), (3, 2)])
+  @test _evaluate_symbolic_local_matrix(
+    duplicate_terms, coefficients, qubit_op_cache_vec[1]
+  ) ≈ 2 * coefficients[2] * qubit_op_cache_vec[1][2].matrix
+
+  fermion_sites = siteinds("Fermion", 1)
+  fermion_op_cache_vec = to_OpCacheVec(fermion_sites, [["I", "C", "Cdag", "N"]])
+  plain_matrix = _evaluate_symbolic_local_matrix(
+    SymbolicLocalMatrix{Int}([(3, 2)]), coefficients, fermion_op_cache_vec[1]
+  )
+  jw_matrix = _evaluate_symbolic_local_matrix(
+    SymbolicLocalMatrix{Int}([(3, -2)]), coefficients, fermion_op_cache_vec[1]
+  )
+
+  expected_plain_matrix = coefficients[2] * fermion_op_cache_vec[1][2].matrix
+  expected_jw_matrix = coefficients[2] * copy(fermion_op_cache_vec[1][2].matrix)
+  expected_jw_matrix[:, 2] .*= -1
+
+  @test plain_matrix ≈ expected_plain_matrix
+  @test jw_matrix ≈ expected_jw_matrix
+
+  return nothing
+end
+
+@testset "SymbolicLocalMatrix" begin
+  test_symbolic_local_matrix_operations()
+end

From 340856c9d30187842c2aa4114a24bca661d820b2 Mon Sep 17 00:00:00 2001
From: Ben Corbett <corbett5@llnl.gov>
Date: Sun, 26 Apr 2026 22:47:26 -0600
Subject: [PATCH 03/16] Step 2

---
 plan.md                             | 24 ++++++++--------
 src/SymbolicLocalMatrix.jl          | 25 +++++++++++++++--
 test/symbolic-local-matrix-tests.jl | 43 ++++++++++++++++++++++++++++-
 3 files changed, 76 insertions(+), 16 deletions(-)

diff --git a/plan.md b/plan.md
index b4ce6f3..afaa0dd 100644
--- a/plan.md
+++ b/plan.md
@@ -34,8 +34,8 @@ H2 = instantiate_MPO(H, sym, coefficients; checkflux=false)
 - Convert each user label `k` to internal symbolic id `k + 1`, reserving
   internal id `1` for the constant one introduced by the vertex-cover algorithm.
 - Carry signs by negating the internal symbolic id. Negative ids can arise from
-  fermionic sorting signs or sign-only basis rewrites; they are not user-facing
-  labels.
+  sign-only basis rewrites and other internal sign handling; they are not
+  user-facing labels.
 - During numeric substitution:
   - internal id `+1` evaluates to `+1.0`,
   - internal id `-1` evaluates to `-1.0`,
@@ -156,22 +156,20 @@ internalize_symbolic_ids!(os::OpIDSum{N,C,Ti}) where {N,C<:Integer,Ti}
 - This helper may mutate the provided `OpIDSum`, matching the existing MPO
   construction behavior where preprocessing, sorting, and duplicate compaction
   are destructive.
-- Map each stored scalar label to a signed internal id in place.
-- Reject scalar label `0`.
+- Map each stored positive scalar label to an internal id in place.
+- Reject scalar labels less than or equal to `0`.
 - Map stored scalar label `k > 0` to internal id `k + 1`.
-- Map stored scalar label `k < 0` to internal id `-(abs(k) + 1)`, treating it
-  as the sign-carrying form of user label `abs(k)`.
 - Keep the same `op_cache_vec`, `abs_tol`, and `modify!`.
 - Add tests proving the in-place conversion maps labels correctly and rejects
-  zero labels.
+  nonpositive labels.
 
 Tests for this step:
 
-- Build a small `OpIDSum{2,Int,Int}` with scalar labels `1`, `2`, and `-3`;
+- Build a small `OpIDSum{2,Int,Int}` with scalar labels `1`, `2`, and `3`;
   after `internalize_symbolic_ids!`, assert that the stored labels are `2`, `3`,
-  and `-4`.
-- Add a term with label `0` and assert `internalize_symbolic_ids!` throws an
-  `ArgumentError`.
+  and `4`.
+- Add terms with labels `0` and `-1` and assert `internalize_symbolic_ids!`
+  throws an `ArgumentError`.
 - Assert the helper is intentionally mutating by checking that the same `os`
   object has updated scalar storage after conversion.
 - Assert `op_cache_vec`, `abs_tol`, and `modify!` are still the original objects
@@ -377,8 +375,8 @@ failure points directly to the layer that broke:
   - unsupported rewrite-factor error.
 - `"internal symbolic ids"`:
   - in-place user-label to internal-id conversion,
-  - zero label rejection,
-  - sign-preserving conversion for labels made negative by preprocessing.
+  - nonpositive label rejection,
+  - negative internal ids from later sign-carrying operations remain supported.
 - `"basis rewrite"`:
   - `+1` rewrite factor preserves id sign,
   - `-1` rewrite factor flips id sign,
diff --git a/src/SymbolicLocalMatrix.jl b/src/SymbolicLocalMatrix.jl
index aecf146..f9d9995 100644
--- a/src/SymbolicLocalMatrix.jl
+++ b/src/SymbolicLocalMatrix.jl
@@ -38,8 +38,29 @@ function _check_symbolic_local_op_id(signed_local_op_id::Integer)::Nothing
 end
 
 function _internal_symbolic_id(user_label::Integer)::Int
-  user_label == 0 && throw(ArgumentError("Symbolic coefficient labels must be nonzero."))
-  return user_label > 0 ? user_label + 1 : -(abs(user_label) + 1)
+  user_label > 0 ||
+    throw(ArgumentError("Symbolic coefficient labels must be greater than zero."))
+  return user_label + 1
+end
+
+"""
+    internalize_symbolic_ids!(os::OpIDSum{N,C,Ti}) where {N,C<:Integer,Ti}
+
+Map the integer scalar labels stored in `os` into the internal symbolic id
+space used by symbolic MPO construction.
+
+User label `k > 0` is stored as `k + 1`, reserving internal id `1` for the
+constant one. Labels less than or equal to zero are invalid. The conversion
+mutates `os` in place and does not alter its operator cache, tolerance, or
+modification callback.
+"""
+function internalize_symbolic_ids!(
+  os::OpIDSum{N,C,Ti}
+)::OpIDSum{N,C,Ti} where {N,C<:Integer,Ti}
+  for i in eachindex(os)
+    os.scalars[i] = C(_internal_symbolic_id(os.scalars[i]))
+  end
+  return os
 end
 
 function _max_user_label(signed_weight_id::Integer)::Int
diff --git a/test/symbolic-local-matrix-tests.jl b/test/symbolic-local-matrix-tests.jl
index 56d76dd..f65342e 100644
--- a/test/symbolic-local-matrix-tests.jl
+++ b/test/symbolic-local-matrix-tests.jl
@@ -7,7 +7,8 @@ using ITensorMPOConstruction:
   _max_user_label,
   _normalize_symbolic_local_matrix!,
   _scale_symbolic_weight,
-  _weight_value
+  _weight_value,
+  internalize_symbolic_ids!
 using ITensors
 using ITensorMPS
 using Test
@@ -24,6 +25,7 @@ function test_symbolic_local_matrix_operations()::Nothing
   @test _internal_symbolic_id(1) == 2
   @test _internal_symbolic_id(3) == 4
   @test_throws ArgumentError _internal_symbolic_id(0)
+  @test_throws ArgumentError _internal_symbolic_id(-1)
 
   @test _scale_symbolic_weight(3, 1) == 3
   @test _scale_symbolic_weight(3, -1) == -3
@@ -74,6 +76,45 @@ function test_symbolic_local_matrix_operations()::Nothing
   return nothing
 end
 
+function test_internalize_symbolic_ids()::Nothing
+  sites = siteinds("Qubit", 2)
+  op_cache_vec = to_OpCacheVec(sites, [["I", "X"], ["I", "Z"]])
+
+  X(n) = OpID(2, n)
+  Z(n) = OpID(2, n)
+
+  modify!(ops) = 1
+  os = OpIDSum{2,Int,Int}(3, op_cache_vec, modify!; abs_tol=0.25)
+  add!(os, 1, X(1))
+  add!(os, 2, Z(2))
+  add!(os, 3, X(1), Z(2))
+
+  original_op_cache_vec = os.op_cache_vec
+  original_abs_tol = os.abs_tol
+  original_modify! = os.modify!
+  original_scalars = os.scalars
+
+  @test internalize_symbolic_ids!(os) === os
+  @test original_scalars === os.scalars
+  @test os.scalars[1:3] == [2, 3, 4]
+  @test os.op_cache_vec === original_op_cache_vec
+  @test os.abs_tol == original_abs_tol
+  @test os.modify! === original_modify!
+
+  zero_label_os = OpIDSum{2,Int,Int}(1, op_cache_vec)
+  add!(zero_label_os, 1, X(1))
+  zero_label_os.scalars[1] = 0
+  @test_throws ArgumentError internalize_symbolic_ids!(zero_label_os)
+
+  negative_label_os = OpIDSum{2,Int,Int}(1, op_cache_vec)
+  add!(negative_label_os, -1, X(1))
+  @test_throws ArgumentError internalize_symbolic_ids!(negative_label_os)
+
+  return nothing
+end
+
 @testset "SymbolicLocalMatrix" begin
   test_symbolic_local_matrix_operations()
+
+  test_internalize_symbolic_ids()
 end

From f97f50a9a669a5355e747439af2fb4efe8005c60 Mon Sep 17 00:00:00 2001
From: Ben Corbett <corbett5@llnl.gov>
Date: Sun, 26 Apr 2026 23:08:00 -0600
Subject: [PATCH 04/16] Step 3

---
 plan.md                             | 12 +++--
 src/OpIDSum.jl                      | 45 ++++++++++++----
 test/symbolic-local-matrix-tests.jl | 82 ++++++++++++++++++++++++++++-
 3 files changed, 125 insertions(+), 14 deletions(-)

diff --git a/plan.md b/plan.md
index afaa0dd..e550da2 100644
--- a/plan.md
+++ b/plan.md
@@ -177,12 +177,16 @@ Tests for this step:
 
 ### 3. Make basis rewrite compatible with internal signed ids
 
-- Reuse `prepare_opID_sum!(os, basis_op_cache_vec)` after the in-place
-  internal-id conversion.
+- Reuse `prepare_opID_sum!(os, basis_op_cache_vec; symbolic_coefficients=true)`
+  after the in-place internal-id conversion. The keyword keeps ordinary numeric
+  `MPO_new` preprocessing on the existing coefficient-scaling path while the
+  symbolic path treats integer coefficients as internal signed ids.
 - Ensure `rewrite_in_operator_basis!` can multiply internal signed ids by basis
   rewrite factors only when the factor is exactly `+1` or `-1`.
-- Reject `im`, `2im`, `0.5`, `2`, and other unsupported factors with a clear
-  `ArgumentError`.
+- Reject `im`, `2im`, `0.5`, `2`, and other unsupported factors with an
+  `ArgumentError`; this may be reported through Julia's threaded task failure
+  wrapper because symbolic rewrites use the same threaded loop as numeric
+  rewrites.
 - Add a sign-only rewrite regression where a same-site product produces `-1`
   and verifies the internal id sign flips while the user label is preserved.
 
diff --git a/src/OpIDSum.jl b/src/OpIDSum.jl
index fc4c36e..57a289e 100644
--- a/src/OpIDSum.jl
+++ b/src/OpIDSum.jl
@@ -384,7 +384,9 @@ function for_equal_sites(f::Function, ops::AbstractVector{<:OpID})::Nothing
 end
 
 """
-    rewrite_in_operator_basis!(os::OpIDSum, basis_op_cache_vec::OpCacheVec)
+    rewrite_in_operator_basis!(
+      os::OpIDSum, basis_op_cache_vec::OpCacheVec; symbolic_coefficients=false
+    ) -> Nothing
 
 Rewrite each same-site operator product in `os` into a single operator drawn
 from `basis_op_cache_vec`.
@@ -397,10 +399,23 @@ matrix, the whole term is set to zero. An error is thrown if a nonzero product
 cannot be represented by a single basis operator. After rewriting, each
 nonzero term is resorted by site and `os.op_cache_vec` is replaced by
 `basis_op_cache_vec`.
+
+When `symbolic_coefficients=true`, the stored coefficients are interpreted as
+internal signed symbolic ids. In that mode, rewrite factors must be exactly
+`+1` or `-1`; unsupported factors throw an `ArgumentError` instead of
+converting the symbolic id into an ordinary numeric coefficient.
 """
 @timeit function rewrite_in_operator_basis!(
-  os::OpIDSum{N,C,Ti}, basis_op_cache_vec::OpCacheVec
-) where {N,C,Ti}
+  os::OpIDSum{N,C,Ti}, basis_op_cache_vec::OpCacheVec; symbolic_coefficients::Bool=false
+)::Nothing where {N,C,Ti}
+  if symbolic_coefficients && !(C <: Integer)
+    throw(
+      ArgumentError(
+        "Symbolic operator-basis rewrites require integer internal coefficient ids.",
+      ),
+    )
+  end
+
   op_cache_vec = os.op_cache_vec
 
   function scale_by_first_nz!(matrix::Matrix)
@@ -456,6 +471,11 @@ nonzero term is resorted by site and `os.op_cache_vec` is replaced by
     end
   end
 
+  function apply_rewrite_factor(scalar, coeff)
+    symbolic_coefficients && return C(_scale_symbolic_weight(scalar, coeff))
+    return scalar * coeff
+  end
+
   Threads.@threads for i in eachindex(os)
     scalar, ops = os[i]
     term_is_zero = false
@@ -484,7 +504,7 @@ nonzero term is resorted by site and `os.op_cache_vec` is replaced by
         )
       end
 
-      scalar *= coeff
+      scalar = apply_rewrite_factor(scalar, coeff)
 
       ops[a] = OpID{Ti}(basis_id, ops[a].n)
       for k in (a + 1):b
@@ -498,7 +518,7 @@ nonzero term is resorted by site and `os.op_cache_vec` is replaced by
 
   os.op_cache_vec = basis_op_cache_vec
 
-  return os
+  return nothing
 end
 
 """
@@ -620,18 +640,25 @@ term has even fermion parity.
 end
 
 """
-    prepare_opID_sum!(os::OpIDSum, basis_op_cache_vec::Union{Nothing,OpCacheVec}) -> Nothing
+    prepare_opID_sum!(
+      os::OpIDSum, basis_op_cache_vec::Union{Nothing,OpCacheVec};
+      symbolic_coefficients=false
+    ) -> Nothing
 
 Apply optional preprocessing to `os` before MPO construction.
 
 Currently this rewrites terms into `basis_op_cache_vec` when a basis cache is
-provided and otherwise leaves `os` unchanged.
+provided and otherwise leaves `os` unchanged. Set `symbolic_coefficients=true`
+when integer coefficients are internal signed symbolic ids rather than numeric
+scalars.
 """
 function prepare_opID_sum!(
-  os::OpIDSum, basis_op_cache_vec::Union{Nothing,OpCacheVec}
+  os::OpIDSum,
+  basis_op_cache_vec::Union{Nothing,OpCacheVec};
+  symbolic_coefficients::Bool=false,
 )::Nothing
   if !isnothing(basis_op_cache_vec)
-    rewrite_in_operator_basis!(os, basis_op_cache_vec)
+    rewrite_in_operator_basis!(os, basis_op_cache_vec; symbolic_coefficients)
   end
 
   return nothing
diff --git a/test/symbolic-local-matrix-tests.jl b/test/symbolic-local-matrix-tests.jl
index f65342e..1da34b5 100644
--- a/test/symbolic-local-matrix-tests.jl
+++ b/test/symbolic-local-matrix-tests.jl
@@ -8,7 +8,8 @@ using ITensorMPOConstruction:
   _normalize_symbolic_local_matrix!,
   _scale_symbolic_weight,
   _weight_value,
-  internalize_symbolic_ids!
+  internalize_symbolic_ids!,
+  prepare_opID_sum!
 using ITensors
 using ITensorMPS
 using Test
@@ -113,8 +114,87 @@ function test_internalize_symbolic_ids()::Nothing
   return nothing
 end
 
+function symbolic_basis_rewrite_cache_vecs(factor)
+  identity_matrix = [1.0 0.0; 0.0 1.0]
+  x_matrix = [0.0 1.0; 1.0 0.0]
+
+  op_cache_vec = [[
+    OpInfo("I", identity_matrix, false, QN()),
+    OpInfo("scaled I", factor * identity_matrix, false, QN()),
+    OpInfo("X", x_matrix, false, QN()),
+  ]]
+  basis_op_cache_vec = [[
+    OpInfo("I", identity_matrix, false, QN()), OpInfo("X", x_matrix, false, QN())
+  ]]
+
+  return op_cache_vec, basis_op_cache_vec
+end
+
+function symbolic_rewrite_opID_sum(factor, user_label::Int=7)
+  op_cache_vec, basis_op_cache_vec = symbolic_basis_rewrite_cache_vecs(factor)
+
+  scaled_identity(n) = OpID(2, n)
+  X(n) = OpID(3, n)
+
+  os = OpIDSum{2,Int,Int}(1, op_cache_vec)
+  add!(os, user_label, scaled_identity(1), X(1))
+  internalize_symbolic_ids!(os)
+
+  return os, basis_op_cache_vec
+end
+
+function test_symbolic_basis_rewrite_sign_factors()::Nothing
+  user_label = 7
+  internal_id = user_label + 1
+
+  os, basis_op_cache_vec = symbolic_rewrite_opID_sum(-1.0, user_label)
+  prepare_opID_sum!(os, basis_op_cache_vec; symbolic_coefficients=true)
+
+  scalar, ops = os[1]
+  nonzero_ops = [op for op in ops if op != zero(op)]
+  @test scalar == -internal_id
+  @test nonzero_ops == [OpID(2, 1)]
+  @test count(op -> op == zero(op), ops) == 1
+  @test os.op_cache_vec === basis_op_cache_vec
+
+  os, basis_op_cache_vec = symbolic_rewrite_opID_sum(1.0, user_label)
+  prepare_opID_sum!(os, basis_op_cache_vec; symbolic_coefficients=true)
+
+  scalar, ops = os[1]
+  nonzero_ops = [op for op in ops if op != zero(op)]
+  @test scalar == internal_id
+  @test nonzero_ops == [OpID(2, 1)]
+  @test count(op -> op == zero(op), ops) == 1
+  @test os.op_cache_vec === basis_op_cache_vec
+
+  return nothing
+end
+
+function test_symbolic_basis_rewrite_rejects_unsupported_factors()::Nothing
+  for factor in (im, 2im, 0.5, 2.0)
+    os, basis_op_cache_vec = symbolic_rewrite_opID_sum(factor)
+    err = nothing
+    try
+      prepare_opID_sum!(os, basis_op_cache_vec; symbolic_coefficients=true)
+    catch e
+      err = e
+    end
+
+    @test !isnothing(err)
+    error_message = sprint(showerror, err)
+    @test occursin("ArgumentError", error_message)
+    @test occursin("Unsupported symbolic rewrite factor", error_message)
+  end
+
+  return nothing
+end
+
 @testset "SymbolicLocalMatrix" begin
   test_symbolic_local_matrix_operations()
 
   test_internalize_symbolic_ids()
+
+  test_symbolic_basis_rewrite_sign_factors()
+
+  test_symbolic_basis_rewrite_rejects_unsupported_factors()
 end

From d53a29576128367c0eb6d6f42c89acc5746350d4 Mon Sep 17 00:00:00 2001
From: Ben Corbett <corbett5@llnl.gov>
Date: Mon, 27 Apr 2026 11:10:06 -0600
Subject: [PATCH 05/16] Step 4

---
 plan.md                             |  18 ++--
 src/large-graph-mpo-common.jl       |  42 +++++---
 test/symbolic-local-matrix-tests.jl | 151 ++++++++++++++++++++++++++++
 3 files changed, 191 insertions(+), 20 deletions(-)

diff --git a/plan.md b/plan.md
index e550da2..3ebd61a 100644
--- a/plan.md
+++ b/plan.md
@@ -220,16 +220,15 @@ Tests for this step:
 Tests for this step:
 
 - Build an `OpIDSum{2,Int,Int}` with two identical operator tuples labeled `1`
-  and `2`; after symbolic graph/block construction, assert the resulting
-  symbolic local matrix contains both signed ids, not their numeric sum.
+  and `2`; after symbolic graph construction and duplicate-right-vertex
+  compaction, assert the resulting graph preserves both signed ids as separate
+  edge weights, not their numeric sum.
 - Add two identical operator tuples with the same label `1`; instantiate with
-  `coefficients[1] = a` and verify the numeric MPO contribution is `2a`.
+  `coefficients[1] = a` through the symbolic local matrix helper and verify the
+  local contribution is `2a`.
 - Add two identical operator tuples whose labels become opposite signed internal
-  ids after preprocessing; verify the symbolic local matrix cancels or
-  instantiates to zero for that contribution.
-- Run the same duplicate-label test through `MPO_symbolic` and compare the
-  instantiated MPO against numeric `MPO_new` built from the corresponding
-  numeric coefficients.
+  ids after preprocessing; verify the derived symbolic local matrix cancels or
+  evaluates to zero for that contribution.
 
 ### 5. Build `SymbolicMPO`
 
@@ -309,6 +308,9 @@ Tests for this step:
 - Instantiate a small symbolic MPO with two different coefficient vectors and
   compare each result against numeric `MPO_new` built from an equivalent numeric
   `OpIDSum`.
+- Run the duplicate-label symbolic graph regression through `MPO_symbolic` and
+  compare the instantiated MPO against numeric `MPO_new` built from the
+  corresponding numeric coefficients.
 - Instantiate with `splitblocks=false` and `splitblocks=true` for a QN-conserving
   fixture and compare both against numeric `MPO_new`.
 - Instantiate with a coefficient vector containing zeros and assert the link
diff --git a/src/large-graph-mpo-common.jl b/src/large-graph-mpo-common.jl
index 839c79c..8732403 100644
--- a/src/large-graph-mpo-common.jl
+++ b/src/large-graph-mpo-common.jl
@@ -200,7 +200,7 @@ function add_to_next_graph!(
 end
 
 """
-    MPOGraph(os::OpIDSum{N,C,Ti}) -> MPOGraph{N,C,Ti}
+    MPOGraph(os::OpIDSum{N,C,Ti}; symbolic_coefficients=false) -> MPOGraph{N,C,Ti}
 
 Convert an `OpIDSum` into the initial bipartite graph used by the
 MPO construction algorithm.
@@ -208,13 +208,21 @@ MPO construction algorithm.
 Operators within each term are put in reverse order (by decreasing site), then
 the terms are sorted along with the scalars. This sorting puts terms that share
 a terminating sequence of operators (which is now at the front of the term) nearby.
-Duplicate terms are then combined, and resulting terms with a weight below
-`os.abs_tol` are dropped. The returned graph is split before the first site:
-left vertices represent the identity incoming link and the operator emitted at
-site 1, while right vertices carry the remaining operator tuples.
+For numeric construction, duplicate terms are then combined and resulting terms
+with a weight below `os.abs_tol` are dropped. When `symbolic_coefficients=true`,
+scalars are interpreted as signed internal symbolic ids, so duplicate operator
+tuples are preserved as separate edge entries instead of being added together.
+The returned graph is split before the first site: left vertices represent the
+identity incoming link and the operator emitted at site 1, while right vertices
+carry the remaining operator tuples.
 """
-@timeit function MPOGraph(os::OpIDSum{N,C,Ti})::MPOGraph{N,C,Ti} where {N,C,Ti}
+@timeit function MPOGraph(
+  os::OpIDSum{N,C,Ti}; symbolic_coefficients::Bool=false
+)::MPOGraph{N,C,Ti} where {N,C,Ti}
   @assert size(os.terms, 1) == N
+  if symbolic_coefficients && !(C <: Integer)
+    throw(ArgumentError("Symbolic MPO graph construction requires integer coefficient ids."))
+  end
 
   ## Reverse the terms in the sum, ignoring trailing identity operators.
   Threads.@threads for i in 1:length(os)
@@ -234,17 +242,27 @@ site 1, while right vertices carry the remaining operator tuples.
     alg=(Threads.nthreads() > 1) ? ThreadsX.QuickSort : Base.QuickSort,
   )
 
-  ## Combine duplicate terms and remove terms which are below the tolerance.
+  ## Combine duplicate numeric terms and remove terms which are below the tolerance.
+  ## Symbolic ids are labels, not numeric weights, so preserve duplicate entries.
   nnz = 0
-  for i in eachindex(os)
-    if i < length(os) && os._data[i] == os._data[i + 1]
-      os.scalars[i + 1] += os.scalars[i]
-      os.scalars[i] = 0
-    elseif abs(os.scalars[i]) > os.abs_tol
+  if symbolic_coefficients
+    for i in eachindex(os)
+      _check_symbolic_weight_id(os.scalars[i])
       nnz += 1
       os.scalars[nnz] = os.scalars[i]
       os._data[nnz] = os._data[i]
     end
+  else
+    for i in eachindex(os)
+      if i < length(os) && os._data[i] == os._data[i + 1]
+        os.scalars[i + 1] += os.scalars[i]
+        os.scalars[i] = 0
+      elseif abs(os.scalars[i]) > os.abs_tol
+        nnz += 1
+        os.scalars[nnz] = os.scalars[i]
+        os._data[nnz] = os._data[i]
+      end
+    end
   end
 
   os.num_terms[] = nnz
diff --git a/test/symbolic-local-matrix-tests.jl b/test/symbolic-local-matrix-tests.jl
index 1da34b5..da702e0 100644
--- a/test/symbolic-local-matrix-tests.jl
+++ b/test/symbolic-local-matrix-tests.jl
@@ -1,5 +1,6 @@
 using ITensorMPOConstruction
 using ITensorMPOConstruction:
+  MPOGraph,
   SymbolicLocalMatrix,
   _append_symbolic_local_matrix_term!,
   _evaluate_symbolic_local_matrix,
@@ -8,8 +9,12 @@ using ITensorMPOConstruction:
   _normalize_symbolic_local_matrix!,
   _scale_symbolic_weight,
   _weight_value,
+  combine_duplicate_adjacent_right_vertices!,
   internalize_symbolic_ids!,
+  left_size,
+  left_vertex,
   prepare_opID_sum!
+using ITensorMPOConstruction: right_size, terms_eq_from
 using ITensors
 using ITensorMPS
 using Test
@@ -189,6 +194,144 @@ function test_symbolic_basis_rewrite_rejects_unsupported_factors()::Nothing
   return nothing
 end
 
+function two_site_qubit_symbolic_fixture()
+  sites = siteinds("Qubit", 2)
+  op_cache_vec = to_OpCacheVec(sites, [["I", "X"] for _ in eachindex(sites)])
+  return sites, op_cache_vec
+end
+
+function symbolic_terms_from_first_left_vertex(g)::SymbolicLocalMatrix{Int}
+  @test left_size(g) == 1
+  lv = left_vertex(g, 1)
+  signed_local_op_id = lv.needs_JW_string ? -lv.op_id : lv.op_id
+  terms = SymbolicLocalMatrix{Int}()
+  for signed_weight_id in g.edge_weights_from_left[1]
+    _append_symbolic_local_matrix_term!(terms, signed_weight_id, signed_local_op_id)
+  end
+  return _normalize_symbolic_local_matrix!(terms)
+end
+
+function compact_duplicate_symbolic_right_vertices!(g)::Nothing
+  combine_duplicate_adjacent_right_vertices!(g, terms_eq_from(2))
+  @test right_size(g) == 1
+  @test g.right_vertex_ids_from_left == [[1 for _ in g.edge_weights_from_left[1]]]
+  return nothing
+end
+
+function test_numeric_mpo_graph_still_compacts_duplicate_terms()::Nothing
+  sites = siteinds("Qubit", 2)
+  op_cache_vec = to_OpCacheVec(sites, [["I", "X", "Z"] for _ in eachindex(sites)])
+
+  X(n) = OpID(2, n)
+  Z(n) = OpID(3, n)
+
+  os = OpIDSum{2,Float64,Int}(4, op_cache_vec; abs_tol=0.25)
+  add!(os, 1.0, X(1), X(2))
+  add!(os, 2.0, X(1), X(2))
+  add!(os, 0.5, Z(1), Z(2))
+  add!(os, -0.4, Z(1), Z(2))
+
+  g = MPOGraph(os)
+
+  @test left_size(g) == 1
+  @test right_size(g) == 1
+  @test left_vertex(g, 1).op_id == 2
+  @test g.right_vertex_ids_from_left == [[1]]
+  @test g.edge_weights_from_left == [[3.0]]
+
+  return nothing
+end
+
+function test_symbolic_mpo_graph_preserves_duplicate_signed_ids()::Nothing
+  _, op_cache_vec = two_site_qubit_symbolic_fixture()
+
+  X(n) = OpID(2, n)
+
+  os = OpIDSum{2,Int,Int}(2, op_cache_vec)
+  add!(os, 1, X(1), X(2))
+  add!(os, 2, X(1), X(2))
+  internalize_symbolic_ids!(os)
+
+  g = MPOGraph(os; symbolic_coefficients=true)
+  @test length(g.edge_weights_from_left[1]) == 2
+  @test sort(g.edge_weights_from_left[1]) == [2, 3]
+
+  compact_duplicate_symbolic_right_vertices!(g)
+  terms = symbolic_terms_from_first_left_vertex(g)
+  @test terms == [(2, 2), (3, 2)]
+
+  coefficients = [11.0, 17.0]
+  evaluated = _evaluate_symbolic_local_matrix(terms, coefficients, op_cache_vec[1])
+  @test evaluated ≈ sum(coefficients) * op_cache_vec[1][2].matrix
+
+  return nothing
+end
+
+function test_symbolic_mpo_graph_preserves_duplicate_label_multiplicity()::Nothing
+  _, op_cache_vec = two_site_qubit_symbolic_fixture()
+
+  X(n) = OpID(2, n)
+
+  os = OpIDSum{2,Int,Int}(2, op_cache_vec)
+  add!(os, 1, X(1), X(2))
+  add!(os, 1, X(1), X(2))
+  internalize_symbolic_ids!(os)
+
+  g = MPOGraph(os; symbolic_coefficients=true)
+  compact_duplicate_symbolic_right_vertices!(g)
+  terms = symbolic_terms_from_first_left_vertex(g)
+  @test terms == [(2, 2), (2, 2)]
+
+  coefficient = 13.0
+  evaluated = _evaluate_symbolic_local_matrix(terms, [coefficient], op_cache_vec[1])
+  @test evaluated ≈ 2 * coefficient * op_cache_vec[1][2].matrix
+
+  return nothing
+end
+
+function signed_symbolic_rewrite_cache_vecs()
+  identity_matrix = [1.0 0.0; 0.0 1.0]
+  x_matrix = [0.0 1.0; 1.0 0.0]
+
+  op_cache_vec = [
+    [
+      OpInfo("I", identity_matrix, false, QN()),
+      OpInfo("X", x_matrix, false, QN()),
+      OpInfo("-X", -x_matrix, false, QN()),
+    ],
+    [OpInfo("I", identity_matrix, false, QN()), OpInfo("X", x_matrix, false, QN())],
+  ]
+  basis_op_cache_vec = [
+    [OpInfo("I", identity_matrix, false, QN()), OpInfo("X", x_matrix, false, QN())],
+    [OpInfo("I", identity_matrix, false, QN()), OpInfo("X", x_matrix, false, QN())],
+  ]
+
+  return op_cache_vec, basis_op_cache_vec
+end
+
+function test_symbolic_mpo_graph_cancels_opposite_signed_ids()::Nothing
+  op_cache_vec, basis_op_cache_vec = signed_symbolic_rewrite_cache_vecs()
+
+  X(n) = OpID(2, n)
+  negX(n) = OpID(3, n)
+
+  os = OpIDSum{2,Int,Int}(2, op_cache_vec)
+  add!(os, 1, X(1), X(2))
+  add!(os, 1, negX(1), X(2))
+  internalize_symbolic_ids!(os)
+  prepare_opID_sum!(os, basis_op_cache_vec; symbolic_coefficients=true)
+
+  g = MPOGraph(os; symbolic_coefficients=true)
+  compact_duplicate_symbolic_right_vertices!(g)
+  terms = symbolic_terms_from_first_left_vertex(g)
+  @test isempty(terms)
+
+  evaluated = _evaluate_symbolic_local_matrix(terms, [7.0], basis_op_cache_vec[1])
+  @test iszero(evaluated)
+
+  return nothing
+end
+
 @testset "SymbolicLocalMatrix" begin
   test_symbolic_local_matrix_operations()
 
@@ -197,4 +340,12 @@ end
   test_symbolic_basis_rewrite_sign_factors()
 
   test_symbolic_basis_rewrite_rejects_unsupported_factors()
+
+  test_numeric_mpo_graph_still_compacts_duplicate_terms()
+
+  test_symbolic_mpo_graph_preserves_duplicate_signed_ids()
+
+  test_symbolic_mpo_graph_preserves_duplicate_label_multiplicity()
+
+  test_symbolic_mpo_graph_cancels_opposite_signed_ids()
 end

From 41ed0a838e82fcc962e38f92ecc5c0f4a3f2a314 Mon Sep 17 00:00:00 2001
From: Ben Corbett <corbett5@llnl.gov>
Date: Mon, 27 Apr 2026 11:57:31 -0600
Subject: [PATCH 06/16] Step 5

---
 plan.md                             |   9 ++-
 src/ITensorMPOConstruction.jl       |   3 +-
 src/MPOConstruction.jl              | 110 +++++++++++++++++++++++++++-
 src/SymbolicLocalMatrix.jl          |   8 ++
 src/large-graph-mpo-common.jl       |  17 +++--
 src/large-graph-mpo-vertex-cover.jl | 109 +++++++++++++++++++++++----
 test/symbolic-local-matrix-tests.jl |  85 +++++++++++++++++++++
 7 files changed, 316 insertions(+), 25 deletions(-)

diff --git a/plan.md b/plan.md
index 3ebd61a..4803ff4 100644
--- a/plan.md
+++ b/plan.md
@@ -258,8 +258,13 @@ type instead of using an abstract field.
 - run `prepare_opID_sum!` and optionally `check_os_for_errors`,
 - build the symbolic graph without numerically adding identifier ids,
 - initialize the dummy left link exactly as `MPO_new` does,
-- run a VC-only symbolic construction driver that mirrors the existing
-  vertex-cover path but writes `SymbolicBlockSparseMatrix{Ti}` blocks,
+- call the existing `resume_MPO_construction!` with `alg="VC"` and symbolic
+  block storage, so symbolic construction reuses `at_site!` and
+  `process_vertex_cover!` instead of maintaining parallel symbolic driver
+  functions,
+- use small storage-dispatch helpers inside the vertex-cover path so numeric
+  blocks still accumulate dense matrices while symbolic blocks append
+  `(signed_weight_id, signed_local_op_id)` terms,
 - when adding a local contribution, append
   `(signed_weight_id, signed_local_op_id)` to the symbolic block instead of
   expanding the dense local operator matrix,
diff --git a/src/ITensorMPOConstruction.jl b/src/ITensorMPOConstruction.jl
index b3b37bb..55b97b8 100644
--- a/src/ITensorMPOConstruction.jl
+++ b/src/ITensorMPOConstruction.jl
@@ -37,6 +37,7 @@ include("MPOConstruction.jl")
 export OpInfo, OpCacheVec, to_OpCacheVec, OpID, OpIDSum
 
 # MPOConstruction.jl
-export resume_MPO_construction!, MPO_new, instantiate_MPO, sparsity, block2_nnz
+export resume_MPO_construction!,
+  MPO_new, SymbolicMPO, MPO_symbolic, instantiate_MPO, sparsity, block2_nnz
 
 end
diff --git a/src/MPOConstruction.jl b/src/MPOConstruction.jl
index 5a8edd8..3f5b961 100644
--- a/src/MPOConstruction.jl
+++ b/src/MPOConstruction.jl
@@ -45,7 +45,7 @@ $_resume_MPO_kwargs
 function resume_MPO_construction!(
   n_init::Int,
   offsets::Vector{Vector{Int}},
-  block_sparse_matrices::Vector{Vector{BlockSparseMatrix{ValType}}},
+  block_sparse_matrices::Vector{Vector{MatrixType}},
   sites::Vector{<:Index},
   llinks::Vector{<:Index},
   g::MPOGraph,
@@ -57,14 +57,14 @@ function resume_MPO_construction!(
   call_back::Function=(
     cur_site::Int,
     offsets::Vector{Vector{Int}},
-    block_sparse_matrices::Vector{Vector{BlockSparseMatrix{ValType}}},
+    block_sparse_matrices::Vector{Vector{MatrixType}},
     sites::Vector{<:Index},
     llinks::Vector{<:Index},
     cur_graph::MPOGraph,
     op_cache_vec::OpCacheVec,
   ) -> nothing,
   output_level::Int=0,
-)::Nothing where {ValType<:Number}
+)::Nothing where {MatrixType}
   @assert !ITensors.using_auto_fermion() # TODO: This should be fixed.
   @assert tol >= 0
 
@@ -73,7 +73,7 @@ function resume_MPO_construction!(
       "At site $n/$(length(sites)) the graph takes up $(Base.format_bytes(Base.summarysize(g)))",
     )
     @time_if output_level 1 "at_site!" g, offsets[n], block_sparse_matrices[n], llinks[n + 1] = at_site!(
-      ValType,
+      MatrixType,
       g,
       n,
       sites,
@@ -136,6 +136,108 @@ function instantiate_MPO(
   return H
 end
 
+"""
+    SymbolicMPO
+
+Intermediate symbolic MPO representation produced by [`MPO_symbolic`](@ref).
+
+The stored `block_sparse_matrices` have the same structural layout as numeric
+MPO construction, but each block entry is a symbolic local matrix term list.
+`max_user_label` records the largest positive coefficient label needed for
+later numeric substitution.
+"""
+struct SymbolicMPO{
+  Ti<:Integer,Sites<:AbstractVector{<:Index},Links<:AbstractVector{<:Index}
+}
+  offsets::Vector{Vector{Int}}
+  block_sparse_matrices::Vector{Vector{SymbolicBlockSparseMatrix{Ti}}}
+  sites::Sites
+  llinks::Links
+  op_cache_vec::OpCacheVec
+  max_user_label::Int
+end
+
+function _check_MPO_symbolic_kwargs(kwargs)::Nothing
+  if haskey(kwargs, :alg)
+    throw(
+      ArgumentError(
+        "Symbolic MPO construction is always vertex-cover based and does not accept an alg keyword.",
+      ),
+    )
+  end
+
+  if !isempty(kwargs)
+    throw(ArgumentError("Unsupported keyword(s) for MPO_symbolic: $(join(keys(kwargs), ", "))."))
+  end
+
+  return nothing
+end
+
+"""
+    MPO_symbolic(os::OpIDSum, sites; basis_op_cache_vec=nothing,
+      check_for_errors=true, combine_qn_sectors=false, output_level=0) -> SymbolicMPO
+
+Run the vertex-cover MPO construction once for an integer-labeled `OpIDSum`,
+storing symbolic local matrix terms that can be numerically instantiated later.
+
+The input `OpIDSum` is preprocessed in place, matching `MPO_new`: positive user
+coefficient labels are remapped to internal symbolic ids, an optional basis
+rewrite is applied, and duplicate symbolic graph entries are preserved.
+"""
+function MPO_symbolic(
+  os::OpIDSum{N,C,Ti},
+  sites::Vector{<:Index};
+  basis_op_cache_vec=nothing,
+  check_for_errors::Bool=true,
+  combine_qn_sectors::Bool=false,
+  output_level::Int=0,
+  kwargs...,
+)::SymbolicMPO where {N,C,Ti}
+  _check_MPO_symbolic_kwargs(kwargs)
+
+  C <: Integer || throw(
+    ArgumentError("MPO_symbolic requires an OpIDSum with integer coefficient labels."),
+  )
+
+  internalize_symbolic_ids!(os)
+  prepare_opID_sum!(os, to_OpCacheVec(sites, basis_op_cache_vec); symbolic_coefficients=true)
+  check_for_errors && check_os_for_errors(os)
+  max_user_label = _max_user_label(os)
+
+  label = "Constructing symbolic MPOGraph from $(length(os)) terms"
+  @time_if output_level 0 label g = MPOGraph(os; symbolic_coefficients=true)
+
+  llinks = Vector{Index}(undef, length(sites) + 1)
+  if hasqns(sites)
+    llinks[1] = Index(QN() => 1; tags="Link,l=0", dir=ITensors.Out)
+  else
+    llinks[1] = Index(1; tags="Link,l=0")
+  end
+
+  offsets = Vector{Vector{Int}}(undef, length(sites))
+  block_sparse_matrices = Vector{Vector{SymbolicBlockSparseMatrix{Ti}}}(undef, length(sites))
+
+  @time_if output_level 0 "Constructing symbolic MPO terms" resume_MPO_construction!(
+    1,
+    offsets,
+    block_sparse_matrices,
+    sites,
+    llinks,
+    g,
+    os.op_cache_vec;
+    alg="VC",
+    combine_qn_sectors,
+    output_level,
+  )
+
+  return SymbolicMPO(offsets, block_sparse_matrices, sites, llinks, os.op_cache_vec, max_user_label)
+end
+
+function MPO_symbolic(os::OpSum, sites::Vector{<:Index}; kwargs...)
+  _check_MPO_symbolic_kwargs(kwargs)
+  throw(ArgumentError("MPO_symbolic supports OpIDSum input only; OpSum input is not supported."))
+end
+
 @doc """
     MPO_new(ValType, os::OpIDSum, sites; basis_op_cache_vec=nothing,
       check_for_errors=true, checkflux=true, splitblocks=true, output_level=0,
diff --git a/src/SymbolicLocalMatrix.jl b/src/SymbolicLocalMatrix.jl
index f9d9995..580033b 100644
--- a/src/SymbolicLocalMatrix.jl
+++ b/src/SymbolicLocalMatrix.jl
@@ -76,6 +76,14 @@ function _max_user_label(terms::SymbolicLocalMatrix)::Int
   return max_user_label
 end
 
+function _max_user_label(os::OpIDSum)::Int
+  max_user_label = 0
+  for i in eachindex(os)
+    max_user_label = max(max_user_label, _max_user_label(os.scalars[i]))
+  end
+  return max_user_label
+end
+
 function _append_symbolic_local_matrix_term!(
   terms::SymbolicLocalMatrix{Ti},
   signed_weight_id::Integer,
diff --git a/src/large-graph-mpo-common.jl b/src/large-graph-mpo-common.jl
index 8732403..10b9a4c 100644
--- a/src/large-graph-mpo-common.jl
+++ b/src/large-graph-mpo-common.jl
@@ -349,6 +349,14 @@ function merge_qn_sectors(
   return new_order, new_qi
 end
 
+function _check_qr_block_storage(::Type{SymbolicBlockSparseMatrix{Ti}})::Nothing where {Ti}
+  throw(ArgumentError("QR construction is only supported for numeric block storage."))
+end
+
+function _check_qr_block_storage(::Type)::Nothing
+  return nothing
+end
+
 """
     at_site!(ValType, g, n, sites, tol, absolute_tol, op_cache_vec, alg;
              combine_qn_sectors, output_level=0)
@@ -370,7 +378,7 @@ merges adjacent outgoing link sectors with the same QN after component ranks are
 known.
 """
 @timeit function at_site!(
-  ::Type{ValType},
+  ::Type{MatrixType},
   g::MPOGraph{N,C,Ti},
   n::Int,
   sites::Vector{<:Index},
@@ -380,9 +388,7 @@ known.
   alg::String;
   combine_qn_sectors::Bool,
   output_level::Int=0,
-)::Tuple{
-  MPOGraph{N,C,Ti},Vector{Int},Vector{BlockSparseMatrix{ValType}},Index
-} where {ValType<:Number,N,C,Ti}
+)::Tuple{MPOGraph{N,C,Ti},Vector{Int},Vector{MatrixType},Index} where {MatrixType,N,C,Ti}
   has_qns = hasqns(sites)
 
   workspace = combine_duplicate_adjacent_right_vertices!(g, terms_eq_from(n + 1))
@@ -395,7 +401,7 @@ known.
   rank_of_cc = zeros(Int, nccs)
 
   ## The MPO tensor for each component.
-  matrix_of_cc = [BlockSparseMatrix{ValType}() for _ in 1:nccs]
+  matrix_of_cc = [_vertex_cover_matrix(MatrixType, 0) for _ in 1:nccs]
 
   ## The QN of each component
   qi_of_cc = Vector{Pair{QN,Int}}(undef, nccs)
@@ -416,6 +422,7 @@ known.
   )
 
   if alg == "QR"
+    _check_qr_block_storage(MatrixType)
     process_qr(
       matrix_of_cc,
       rank_of_cc,
diff --git a/src/large-graph-mpo-vertex-cover.jl b/src/large-graph-mpo-vertex-cover.jl
index 3a548af..65b55e2 100644
--- a/src/large-graph-mpo-vertex-cover.jl
+++ b/src/large-graph-mpo-vertex-cover.jl
@@ -1,3 +1,92 @@
+function _vertex_cover_matrix(
+  ::Type{BlockSparseMatrix{ValType}}, rank::Int
+)::BlockSparseMatrix{ValType} where {ValType}
+  return [Dictionary{Int,Matrix{ValType}}() for _ in 1:rank]
+end
+
+function _vertex_cover_matrix(
+  ::Type{SymbolicBlockSparseMatrix{Ti}}, rank::Int
+)::SymbolicBlockSparseMatrix{Ti} where {Ti<:Integer}
+  return [Dictionary{Int,SymbolicLocalMatrix{Ti}}() for _ in 1:rank]
+end
+
+function _signed_symbolic_local_op_id(lv::LeftVertex, ::Type{Ti})::Ti where {Ti<:Integer}
+  signed_local_op_id = lv.needs_JW_string ? -Int(lv.op_id) : Int(lv.op_id)
+  return Ti(signed_local_op_id)
+end
+
+function _add_vertex_cover_term!(
+  matrix::BlockSparseMatrix{ValType},
+  m::Int,
+  lv::LeftVertex,
+  weight::Number,
+  op_cache::Vector{OpInfo},
+  site_dim::Int,
+)::Nothing where {ValType}
+  local_op = op_cache[lv.op_id].matrix
+  matrix_element = get!(matrix[m], lv.link) do
+    zeros(ValType, site_dim, site_dim)
+  end
+  add_to_local_matrix!(matrix_element, weight, local_op, lv.needs_JW_string)
+  return nothing
+end
+
+function _add_vertex_cover_term!(
+  matrix::SymbolicBlockSparseMatrix{Ti},
+  m::Int,
+  lv::LeftVertex,
+  signed_weight_id::Integer,
+  op_cache::Vector{OpInfo},
+  site_dim::Int,
+)::Nothing where {Ti<:Integer}
+  terms = get!(matrix[m], lv.link) do
+    SymbolicLocalMatrix{Ti}()
+  end
+  _append_symbolic_local_matrix_term!(
+    terms, signed_weight_id, _signed_symbolic_local_op_id(lv, Ti)
+  )
+  return nothing
+end
+
+function _add_vertex_cover_unit_term!(
+  matrix::BlockSparseMatrix{ValType},
+  m::Int,
+  lv::LeftVertex,
+  op_cache::Vector{OpInfo},
+  site_dim::Int,
+  ::Type{C},
+)::Nothing where {ValType,C}
+  _add_vertex_cover_term!(matrix, m, lv, one(ValType), op_cache, site_dim)
+  return nothing
+end
+
+function _add_vertex_cover_unit_term!(
+  matrix::SymbolicBlockSparseMatrix{Ti},
+  m::Int,
+  lv::LeftVertex,
+  op_cache::Vector{OpInfo},
+  site_dim::Int,
+  ::Type{C},
+)::Nothing where {Ti<:Integer,C<:Integer}
+  _add_vertex_cover_term!(matrix, m, lv, one(C), op_cache, site_dim)
+  return nothing
+end
+
+function _finalize_vertex_cover_matrix!(matrix::BlockSparseMatrix)::Nothing
+  return nothing
+end
+
+function _finalize_vertex_cover_matrix!(
+  matrix::SymbolicBlockSparseMatrix
+)::Nothing
+  for block in matrix
+    for terms in values(block)
+      _normalize_symbolic_local_matrix!(terms)
+    end
+  end
+  return nothing
+end
+
 """
     process_vertex_cover!(
       matrix_of_cc, rank_of_cc, next_edges_of_cc, g, ccs, n, sites, op_cache_vec
@@ -12,7 +101,7 @@ stored edge weights. On the final site, only the local tensor blocks are needed;
 otherwise this also builds `next_edges_of_cc` for the graph at site `n + 1`.
 """
 @timeit function process_vertex_cover!(
-  matrix_of_cc::Vector{BlockSparseMatrix{ValType}},
+  matrix_of_cc::Vector{MatrixType},
   rank_of_cc::Vector{Int},
   next_edges_of_cc::Vector{Matrix{Tuple{Vector{Int},Vector{C}}}},
   g::MPOGraph{N,C,Ti},
@@ -20,7 +109,7 @@ otherwise this also builds `next_edges_of_cc` for the graph at site `n + 1`.
   n::Int,
   sites::Vector{<:Index},
   op_cache_vec::OpCacheVec,
-)::Nothing where {ValType<:Number,N,C,Ti}
+)::Nothing where {MatrixType,N,C,Ti}
   site_dim = dim(sites[n])
   op_cache = op_cache_vec[n]
 
@@ -44,18 +133,15 @@ otherwise this also builds `next_edges_of_cc` for the graph at site `n + 1`.
     rank = length(left_cover) + length(right_cover)
     rank_of_cc[cc] = rank
 
-    matrix = [Dictionary{Int,Matrix{ValType}}() for _ in 1:rank]
+    matrix = _vertex_cover_matrix(MatrixType, rank)
     matrix_of_cc[cc] = matrix
 
     ## Construct the tensor from the left cover.
     @inbounds for m in eachindex(left_cover)
       lv_id = lvs_of_component[left_cover[m]]
       lv = left_vertex(g, lv_id)
-      local_op = op_cache[lv.op_id].matrix
 
-      matrix_element = zeros(ValType, site_dim, site_dim)
-      add_to_local_matrix!(matrix_element, one(ValType), local_op, lv.needs_JW_string)
-      set!(matrix[m], lv.link, matrix_element)
+      _add_vertex_cover_unit_term!(matrix, m, lv, op_cache, site_dim, C)
     end
 
     ## Construct the tensor from the right cover.
@@ -82,20 +168,17 @@ otherwise this also builds `next_edges_of_cc` for the graph at site `n + 1`.
 
       @inbounds for lv_id in uncovered_left_ids
         lv = left_vertex(g, lv_id)
-        local_op = op_cache[lv.op_id].matrix
 
         for (rv_id, weight) in weighted_edge_iterator(g, lv_id)
           m = right_cover_m[position_of_rvs_in_component[rv_id]]
 
-          matrix_element = get!(matrix[m], lv.link) do
-            zeros(ValType, site_dim, site_dim)
-          end
-
-          add_to_local_matrix!(matrix_element, weight, local_op, lv.needs_JW_string)
+          _add_vertex_cover_term!(matrix, m, lv, weight, op_cache, site_dim)
         end
       end
     end
 
+    _finalize_vertex_cover_matrix!(matrix)
+
     n == length(sites) && continue
 
     ## Preallocate space for next_edges
diff --git a/test/symbolic-local-matrix-tests.jl b/test/symbolic-local-matrix-tests.jl
index da702e0..25fbccb 100644
--- a/test/symbolic-local-matrix-tests.jl
+++ b/test/symbolic-local-matrix-tests.jl
@@ -1,6 +1,8 @@
 using ITensorMPOConstruction
 using ITensorMPOConstruction:
   MPOGraph,
+  MPO_symbolic,
+  SymbolicMPO,
   SymbolicLocalMatrix,
   _append_symbolic_local_matrix_term!,
   _evaluate_symbolic_local_matrix,
@@ -332,6 +334,83 @@ function test_symbolic_mpo_graph_cancels_opposite_signed_ids()::Nothing
   return nothing
 end
 
+function transverse_field_ising_symbolic_fixture()
+  sites = siteinds("Qubit", 2)
+  op_cache_vec = to_OpCacheVec(sites, [["I", "X", "Z"] for _ in eachindex(sites)])
+
+  X(n) = OpID(2, n)
+  Z(n) = OpID(3, n)
+
+  os = OpIDSum{2,Int,Int}(3, op_cache_vec)
+  add!(os, 1, X(1), X(2))
+  add!(os, 2, Z(1))
+  add!(os, 3, Z(2))
+
+  return os, sites
+end
+
+function symbolic_mpo_terms(sym::SymbolicMPO)::Vector{Tuple{Int,Int}}
+  terms = Tuple{Int,Int}[]
+  for site_matrices in sym.block_sparse_matrices
+    for matrix_of_component in site_matrices
+      for block in matrix_of_component
+        for block_terms in values(block)
+          append!(terms, block_terms)
+        end
+      end
+    end
+  end
+  return terms
+end
+
+function test_symbolic_mpo_construction_metadata()::Nothing
+  os, sites = transverse_field_ising_symbolic_fixture()
+
+  sym = MPO_symbolic(os, sites)
+
+  @test sym isa SymbolicMPO
+  @test length(sym.offsets) == length(sites)
+  @test length(sym.block_sparse_matrices) == length(sites)
+  @test length(sym.llinks) == length(sites) + 1
+  @test sym.max_user_label == 3
+  @test sym.op_cache_vec === os.op_cache_vec
+  @test !isempty(symbolic_mpo_terms(sym))
+
+  return nothing
+end
+
+function test_symbolic_mpo_rejects_unsupported_public_inputs()::Nothing
+  os, sites = transverse_field_ising_symbolic_fixture()
+  @test_throws ArgumentError MPO_symbolic(os, sites; alg="QR")
+
+  noninteger_os = OpIDSum{2,Float64,Int}(1, os.op_cache_vec)
+  add!(noninteger_os, 1.0, OpID(2, 1), OpID(2, 2))
+  @test_throws ArgumentError MPO_symbolic(noninteger_os, sites)
+
+  op_sum = OpSum()
+  op_sum += "X", 1
+  @test_throws ArgumentError MPO_symbolic(op_sum, sites)
+
+  return nothing
+end
+
+function test_symbolic_mpo_uses_negative_local_op_id_for_jw_terms()::Nothing
+  sites = siteinds("Fermion", 2)
+  op_cache_vec = to_OpCacheVec(sites, [["I", "C", "Cdag", "N"] for _ in eachindex(sites)])
+
+  C(n) = OpID(2, n)
+  Cdag(n) = OpID(3, n)
+
+  os = OpIDSum{2,Int,Int}(1, op_cache_vec)
+  add!(os, 1, Cdag(1), C(2))
+
+  sym = MPO_symbolic(os, sites)
+  terms = symbolic_mpo_terms(sym)
+  @test any(term -> term[2] < 0, terms)
+
+  return nothing
+end
+
 @testset "SymbolicLocalMatrix" begin
   test_symbolic_local_matrix_operations()
 
@@ -348,4 +427,10 @@ end
   test_symbolic_mpo_graph_preserves_duplicate_label_multiplicity()
 
   test_symbolic_mpo_graph_cancels_opposite_signed_ids()
+
+  test_symbolic_mpo_construction_metadata()
+
+  test_symbolic_mpo_rejects_unsupported_public_inputs()
+
+  test_symbolic_mpo_uses_negative_local_op_id_for_jw_terms()
 end

From 3571dd69c2f4ebfcf4771c4731994dff0cd90780 Mon Sep 17 00:00:00 2001
From: Ben Corbett <corbett5@llnl.gov>
Date: Mon, 27 Apr 2026 12:24:11 -0600
Subject: [PATCH 07/16] Step 6

---
 plan.md                             |  12 +--
 src/MPOConstruction.jl              | 106 +++++++++++++++++++++++
 test/symbolic-local-matrix-tests.jl | 130 ++++++++++++++++++++++++++++
 3 files changed, 242 insertions(+), 6 deletions(-)

diff --git a/plan.md b/plan.md
index 4803ff4..1ee7024 100644
--- a/plan.md
+++ b/plan.md
@@ -296,17 +296,17 @@ instantiate_MPO(sym::SymbolicMPO, coefficients; splitblocks=true, checkflux=fals
 It should:
 
 - validate the coefficient vector length against `sym.max_user_label`,
-- evaluate symbolic block entries into numeric block matrices using
-  `sym.op_cache_vec`,
-- preserve the symbolic structural pattern even when evaluated values are zero,
+- at each site, evaluate `sym.block_sparse_matrices[n]` with `sym.op_cache_vec`
+  into a `Vector{BlockSparseMatrix{C}}`,
+- preserve the symbolic link layout even when evaluated values are zero,
 - call the existing tensor assembly path with the requested `splitblocks` value,
 - contract away the dummy boundary links the same way as the low-level
   `instantiate_MPO(offsets, block_sparse_matrices, sites, llinks; ...)`,
 - optionally run `ITensors.checkflux`.
 
-For QN tensors, structural block discovery must be based on symbolic local
-matrix terms before numeric substitution, so a coefficient vector containing
-zeros does not change the block structure or link layout.
+For QN tensors, structural block discovery uses the existing numeric
+`to_ITensor` path after symbolic block evaluation. A coefficient vector
+containing zeros must not change the stored symbolic MPO link layout.
 
 Tests for this step:
 
diff --git a/src/MPOConstruction.jl b/src/MPOConstruction.jl
index 3f5b961..032bb18 100644
--- a/src/MPOConstruction.jl
+++ b/src/MPOConstruction.jl
@@ -157,6 +157,112 @@ struct SymbolicMPO{
   max_user_label::Int
 end
 
+function _check_symbolic_coefficients(sym::SymbolicMPO, coefficients::AbstractVector)::Nothing
+  length(coefficients) >= sym.max_user_label || throw(
+    ArgumentError(
+      "Missing coefficient for symbolic label $(sym.max_user_label). " *
+      "Received $(length(coefficients)) coefficients.",
+    ),
+  )
+  return nothing
+end
+
+function _symbolic_mpo_eltype(sym::SymbolicMPO, coefficients::AbstractVector)::Type
+  val_type = promote_type(Float64, eltype(coefficients))
+  for op_cache in sym.op_cache_vec
+    for op_info in op_cache
+      val_type = promote_type(val_type, eltype(op_info.matrix))
+    end
+  end
+  return val_type
+end
+
+function _evaluate_symbolic_block_sparse_matrix(
+  symbolic_matrix::SymbolicBlockSparseMatrix{Ti},
+  coefficients::AbstractVector,
+  op_cache::Vector{OpInfo},
+  ::Type{C},
+)::BlockSparseMatrix{C} where {Ti,C}
+  matrix = [Dictionary{Int,Matrix{C}}() for _ in eachindex(symbolic_matrix)]
+
+  for right_link in eachindex(symbolic_matrix)
+    for (left_link, terms) in pairs(symbolic_matrix[right_link])
+      block = convert(Matrix{C}, _evaluate_symbolic_local_matrix(terms, coefficients, op_cache))
+      set!(matrix[right_link], left_link, block)
+    end
+  end
+
+  return matrix
+end
+
+function _evaluate_symbolic_block_sparse_matrices(
+  symbolic_matrices::Vector{SymbolicBlockSparseMatrix{Ti}},
+  coefficients::AbstractVector,
+  op_cache::Vector{OpInfo},
+  ::Type{C},
+)::Vector{BlockSparseMatrix{C}} where {Ti,C}
+  matrices = Vector{BlockSparseMatrix{C}}(undef, length(symbolic_matrices))
+  for cc in eachindex(symbolic_matrices)
+    matrices[cc] = _evaluate_symbolic_block_sparse_matrix(
+      symbolic_matrices[cc], coefficients, op_cache, C
+    )
+  end
+  return matrices
+end
+
+"""
+    instantiate_MPO(sym::SymbolicMPO, coefficients; splitblocks=true, checkflux=false) -> MPO
+
+Evaluate a symbolic MPO with a numeric coefficient vector and return a fresh
+numeric `MPO`.
+
+The coefficient vector is indexed by the original positive labels used to build
+`sym`. Each site's symbolic block storage is first evaluated into ordinary
+`BlockSparseMatrix` storage, then converted to ITensors through the same
+assembly path used by numeric MPO construction.
+"""
+function instantiate_MPO(
+  sym::SymbolicMPO, coefficients::AbstractVector; splitblocks::Bool=true, checkflux::Bool=false
+)::MPO
+  _check_symbolic_coefficients(sym, coefficients)
+  C = _symbolic_mpo_eltype(sym, coefficients)
+
+  H = MPO(sym.sites)
+
+  @timeit "to_ITensor" Threads.@threads for n in eachindex(sym.sites)
+    block_sparse_matrices = _evaluate_symbolic_block_sparse_matrices(
+      sym.block_sparse_matrices[n], coefficients, sym.op_cache_vec[n], C
+    )
+    H[n] = to_ITensor(
+      sym.offsets[n],
+      block_sparse_matrices,
+      sym.llinks[n],
+      sym.llinks[n + 1],
+      sym.sites[n],
+      ;
+      splitblocks,
+    )
+  end
+
+  # Remove the dummy link indices from the left and right.
+  L = ITensor(sym.llinks[1])
+  L[end] = 1.0
+
+  R = ITensor(dag(sym.llinks[end]))
+  R[1] = 1.0
+
+  H[1] *= L
+  H[length(sym.sites)] *= R
+
+  if checkflux
+    @timeit "checkflux" Threads.@threads for n in eachindex(sym.sites)
+      ITensors.checkflux(H[n])
+    end
+  end
+
+  return H
+end
+
 function _check_MPO_symbolic_kwargs(kwargs)::Nothing
   if haskey(kwargs, :alg)
     throw(
diff --git a/test/symbolic-local-matrix-tests.jl b/test/symbolic-local-matrix-tests.jl
index 25fbccb..6117eee 100644
--- a/test/symbolic-local-matrix-tests.jl
+++ b/test/symbolic-local-matrix-tests.jl
@@ -349,6 +349,60 @@ function transverse_field_ising_symbolic_fixture()
   return os, sites
 end
 
+function transverse_field_ising_numeric_opidsum(sites, coefficients)
+  op_cache_vec = to_OpCacheVec(sites, [["I", "X", "Z"] for _ in eachindex(sites)])
+
+  X(n) = OpID(2, n)
+  Z(n) = OpID(3, n)
+
+  os = OpIDSum{2,eltype(coefficients),Int}(3, op_cache_vec)
+  add!(os, coefficients[1], X(1), X(2))
+  add!(os, coefficients[2], Z(1))
+  add!(os, coefficients[3], Z(2))
+
+  return os
+end
+
+function qn_number_symbolic_fixture()
+  sites = siteinds("Fermion", 2; conserve_qns=true)
+  op_cache_vec = to_OpCacheVec(sites, [["I", "N"] for _ in eachindex(sites)])
+
+  Nop(n) = OpID(2, n)
+
+  os = OpIDSum{2,Int,Int}(3, op_cache_vec)
+  add!(os, 1, Nop(1))
+  add!(os, 2, Nop(2))
+  add!(os, 3, Nop(1), Nop(2))
+
+  return os, sites
+end
+
+function qn_number_numeric_opidsum(sites, coefficients)
+  op_cache_vec = to_OpCacheVec(sites, [["I", "N"] for _ in eachindex(sites)])
+
+  Nop(n) = OpID(2, n)
+
+  os = OpIDSum{2,eltype(coefficients),Int}(3, op_cache_vec)
+  add!(os, coefficients[1], Nop(1))
+  add!(os, coefficients[2], Nop(2))
+  add!(os, coefficients[3], Nop(1), Nop(2))
+
+  return os
+end
+
+function mpo_relative_error(A::MPO, B::MPO)::Float64
+  AmB = add(A, -B; alg="directsum")
+  numerator = real(inner(AmB, AmB))
+  denominator = real(inner(A, A))
+  return iszero(denominator) ? numerator : numerator / denominator
+end
+
+function test_mpos_approx_equal(A::MPO, B::MPO; tol::Real=1e-10)::Nothing
+  relative_error = mpo_relative_error(A, B)
+  @test relative_error < tol
+  return nothing
+end
+
 function symbolic_mpo_terms(sym::SymbolicMPO)::Vector{Tuple{Int,Int}}
   terms = Tuple{Int,Int}[]
   for site_matrices in sym.block_sparse_matrices
@@ -394,6 +448,76 @@ function test_symbolic_mpo_rejects_unsupported_public_inputs()::Nothing
   return nothing
 end
 
+function test_symbolic_mpo_fresh_instantiation_matches_numeric()::Nothing
+  os, sites = transverse_field_ising_symbolic_fixture()
+  sym = MPO_symbolic(os, sites)
+
+  for coefficients in ([1.25, -0.5, 0.75], [-2.0, 0.25, 1.5])
+    symbolic_mpo = instantiate_MPO(sym, coefficients; splitblocks=true)
+    numeric_mpo = MPO_new(
+      transverse_field_ising_numeric_opidsum(sites, coefficients),
+      sites;
+      alg="VC",
+      splitblocks=true,
+      checkflux=false,
+    )
+
+    test_mpos_approx_equal(symbolic_mpo, numeric_mpo)
+  end
+
+  return nothing
+end
+
+function test_symbolic_mpo_duplicate_labels_instantiate_like_numeric()::Nothing
+  sites, op_cache_vec = two_site_qubit_symbolic_fixture()
+
+  X(n) = OpID(2, n)
+
+  symbolic_os = OpIDSum{2,Int,Int}(2, op_cache_vec)
+  add!(symbolic_os, 1, X(1), X(2))
+  add!(symbolic_os, 2, X(1), X(2))
+
+  sym = MPO_symbolic(symbolic_os, sites)
+  coefficients = [2.0, 3.0]
+  symbolic_mpo = instantiate_MPO(sym, coefficients; splitblocks=true)
+
+  numeric_os = OpIDSum{2,Float64,Int}(2, op_cache_vec)
+  add!(numeric_os, coefficients[1], X(1), X(2))
+  add!(numeric_os, coefficients[2], X(1), X(2))
+  numeric_mpo = MPO_new(numeric_os, sites; alg="VC", splitblocks=true, checkflux=false)
+
+  test_mpos_approx_equal(symbolic_mpo, numeric_mpo)
+
+  return nothing
+end
+
+function test_symbolic_mpo_qn_fresh_instantiation()::Nothing
+  os, sites = qn_number_symbolic_fixture()
+  sym = MPO_symbolic(os, sites)
+
+  coefficients = [0.8, -1.1, 0.3]
+  for splitblocks in (false, true)
+    symbolic_mpo = instantiate_MPO(sym, coefficients; splitblocks, checkflux=true)
+    numeric_mpo = MPO_new(
+      qn_number_numeric_opidsum(sites, coefficients),
+      sites;
+      alg="VC",
+      splitblocks,
+      checkflux=true,
+    )
+
+    test_mpos_approx_equal(symbolic_mpo, numeric_mpo)
+  end
+
+  nonzero_mpo = instantiate_MPO(sym, [1.0, 2.0, 3.0]; splitblocks=true)
+  zero_mpo = instantiate_MPO(sym, [0.0, 0.0, 0.0]; splitblocks=true)
+
+  @test linkdims(zero_mpo) == linkdims(nonzero_mpo)
+  @test_throws ArgumentError instantiate_MPO(sym, [1.0, 2.0])
+
+  return nothing
+end
+
 function test_symbolic_mpo_uses_negative_local_op_id_for_jw_terms()::Nothing
   sites = siteinds("Fermion", 2)
   op_cache_vec = to_OpCacheVec(sites, [["I", "C", "Cdag", "N"] for _ in eachindex(sites)])
@@ -432,5 +556,11 @@ end
 
   test_symbolic_mpo_rejects_unsupported_public_inputs()
 
+  test_symbolic_mpo_fresh_instantiation_matches_numeric()
+
+  test_symbolic_mpo_duplicate_labels_instantiate_like_numeric()
+
+  test_symbolic_mpo_qn_fresh_instantiation()
+
   test_symbolic_mpo_uses_negative_local_op_id_for_jw_terms()
 end

From 4d5ee80318203a0e6faba2e35fea517430e076ef Mon Sep 17 00:00:00 2001
From: Ben Corbett <corbett5@llnl.gov>
Date: Mon, 27 Apr 2026 14:37:16 -0600
Subject: [PATCH 08/16] Step 7.

---
 plan.md                             |  18 ++--
 src/ITensorMPOConstruction.jl       |   2 +-
 src/MPOConstruction.jl              | 127 +++++++++++++++++++++++++---
 test/symbolic-local-matrix-tests.jl |  57 ++++++++++++-
 4 files changed, 186 insertions(+), 18 deletions(-)

diff --git a/plan.md b/plan.md
index 1ee7024..59d6e53 100644
--- a/plan.md
+++ b/plan.md
@@ -342,29 +342,37 @@ already fixes the block layout, so evaluation must use the layout of `H` or
 - check `length(H) == length(sym.sites)`,
 - check site indices are compatible with `sym.sites`,
 - check link dimensions match `sym.llinks`,
-- overwrite each tensor in `H` with the evaluated tensor for that site,
-- preserve the existing `MPO` object identity.
+- temporarily add the dummy boundary links back to the first and last MPO
+  tensors by multiplying by one-entry dummy tensors,
+- for each site, zero the existing entries with `H[n] .= 0` and refill through
+  the same four-index coordinate layout used by the numeric tensor assembly,
+- contract away the dummy boundary links again at the end,
+- preserve the existing `MPO` object identity and link layout.
 
 `instantiate_MPO(H_template, sym, coefficients)` should:
 
 - reject incompatible templates with `ArgumentError`,
-- copy or reuse a compatible template,
+- copy a compatible template's tensors before filling them,
 - call `instantiate_MPO!`,
 - return the instantiated MPO.
 
-It is acceptable for V1 to overwrite tensors rather than doing lower-level
-storage mutation, as long as the public behavior and link structure are stable.
+This path should avoid rebuilding ITensors for compatible templates; the
+expected fast path is to reuse the existing structural storage for the nonzero
+entries after temporarily restoring the boundary-link layout.
 
 Tests for this step:
 
 - Build `H = instantiate_MPO(sym, coeffs1; splitblocks=true)`, then call
   `instantiate_MPO!(H, sym, coeffs2)` and compare `H` against
   `instantiate_MPO(sym, coeffs2; splitblocks=true)`.
+- Assert the `MPO` object identity and link dimensions are unchanged after
+  `instantiate_MPO!`.
 - Assert `instantiate_MPO!(H, sym, coeffs; splitblocks=false)` throws a
   `MethodError` or `ArgumentError`, since this overload must not accept
   `splitblocks`.
 - Call `instantiate_MPO(H, sym, coeffs2)` with a compatible template and compare
   against fresh instantiation using the same layout.
+- Assert template-assisted instantiation does not mutate the template MPO.
 - Build an incompatible template with different sites or link dimensions and
   assert `instantiate_MPO(H_bad, sym, coeffs)` throws an `ArgumentError`.
 
diff --git a/src/ITensorMPOConstruction.jl b/src/ITensorMPOConstruction.jl
index 55b97b8..836fa4b 100644
--- a/src/ITensorMPOConstruction.jl
+++ b/src/ITensorMPOConstruction.jl
@@ -38,6 +38,6 @@ export OpInfo, OpCacheVec, to_OpCacheVec, OpID, OpIDSum
 
 # MPOConstruction.jl
 export resume_MPO_construction!,
-  MPO_new, SymbolicMPO, MPO_symbolic, instantiate_MPO, sparsity, block2_nnz
+  MPO_new, SymbolicMPO, MPO_symbolic, instantiate_MPO, instantiate_MPO!, sparsity, block2_nnz
 
 end
diff --git a/src/MPOConstruction.jl b/src/MPOConstruction.jl
index 032bb18..02e6911 100644
--- a/src/MPOConstruction.jl
+++ b/src/MPOConstruction.jl
@@ -91,6 +91,24 @@ function resume_MPO_construction!(
   return nothing
 end
 
+function _remove_symbolic_mpo_boundary_links!(H::MPO, llinks::Vector{<:Index})::Nothing
+  @assert ndims(H[1]) == 4
+  @assert ndims(H[end]) == 4
+
+  L = ITensor(llinks[1])
+  L[end] = 1.0
+  H[1] *= L
+
+  R = ITensor(dag(llinks[end]))
+  R[1] = 1.0
+  H[end] *= R
+
+  @assert ndims(H[1]) == 3
+  @assert ndims(H[end]) == 3
+
+  return nothing
+end
+
 """
     instantiate_MPO(offsets, block_sparse_matrices, sites, llinks; splitblocks, checkflux) -> MPO
 
@@ -117,15 +135,7 @@ function instantiate_MPO(
     )
   end
 
-  # Remove the dummy link indices from the left and right.
-  L = ITensor(llinks[1])
-  L[end] = 1.0
-
-  R = ITensor(dag(llinks[end]))
-  R[1] = 1.0
-
-  H[1] *= L
-  H[length(sites)] *= R
+  _remove_symbolic_mpo_boundary_links!(H, llinks)
 
   if checkflux
     @timeit "checkflux" Threads.@threads for n in 1:length(sites)
@@ -158,9 +168,9 @@ struct SymbolicMPO{
 end
 
 function _check_symbolic_coefficients(sym::SymbolicMPO, coefficients::AbstractVector)::Nothing
-  length(coefficients) >= sym.max_user_label || throw(
+  length(coefficients) == sym.max_user_label || throw(
     ArgumentError(
-      "Missing coefficient for symbolic label $(sym.max_user_label). " *
+      "Expected $(sym.max_user_label) coefficients. " *
       "Received $(length(coefficients)) coefficients.",
     ),
   )
@@ -263,6 +273,101 @@ function instantiate_MPO(
   return H
 end
 
+function _add_symbolic_mpo_boundary_links!(H::MPO, sym::SymbolicMPO)::Nothing
+  L = ITensor(dag(sym.llinks[1]))
+  L[1] = 1.0
+  H[1] *= L
+
+  R = ITensor(sym.llinks[end])
+  R[1] = 1.0
+  H[end] *= R
+
+  return nothing
+end
+
+function _fill_symbolic_mpo_tensor_with_boundary_links!(
+  tensor::ITensor, sym::SymbolicMPO, coefficients::AbstractVector, n::Int
+)::ITensor
+  llink = dag(sym.llinks[n])
+  rlink = sym.llinks[n + 1]
+  site = sym.sites[n]
+  tensor = ITensors.permute(tensor, llink, rlink, prime(site), dag(site); allow_alias=true)
+  tensor .= 0
+
+  op_cache = sym.op_cache_vec[n]
+
+  for (offset, matrix) in zip(sym.offsets[n], sym.block_sparse_matrices[n])
+    for (right_link, column) in enumerate(matrix)
+      shifted_right_link = right_link + offset
+      for (left_link, terms) in pairs(column)
+        local_matrix = _evaluate_symbolic_local_matrix(terms, coefficients, op_cache)
+        for i in axes(local_matrix, 1)
+          for j in axes(local_matrix, 2)
+            iszero(local_matrix[i, j]) && continue
+            tensor[left_link, shifted_right_link, i, j] = local_matrix[i, j]
+          end
+        end
+      end
+    end
+  end
+
+  return tensor
+end
+
+"""
+    instantiate_MPO!(H::MPO, sym::SymbolicMPO, coefficients; checkflux=false) -> MPO
+
+Evaluate `sym` with `coefficients` into the existing MPO `H`.
+
+The target MPO fixes the tensor block layout, so this overload does not accept
+`splitblocks`. Site indices must match `sym`, and internal link dimensions must
+match the symbolic link dimensions. The dummy boundary links are temporarily
+restored so each tensor can be zeroed and refilled through the same four-index
+layout used by numeric tensor assembly, then contracted away again. The MPO
+object `H` is updated in place and returned.
+"""
+function instantiate_MPO!(
+  H::MPO, sym::SymbolicMPO, coefficients::AbstractVector; checkflux::Bool=false
+)::MPO
+  _check_symbolic_coefficients(sym, coefficients)
+
+  if length(H) != length(sym.sites)
+    throw(ArgumentError(
+      "Template MPO length $(length(H)) is incompatible with SymbolicMPO length $(length(sym.sites)).",
+    ))
+  end
+
+  _add_symbolic_mpo_boundary_links!(H, sym)
+
+  Threads.@threads for n in eachindex(sym.sites)
+    H[n] = _fill_symbolic_mpo_tensor_with_boundary_links!(H[n], sym, coefficients, n)
+  end
+
+  _remove_symbolic_mpo_boundary_links!(H, sym.llinks)
+
+  if checkflux
+    @timeit "checkflux" Threads.@threads for n in eachindex(sym.sites)
+      ITensors.checkflux(H[n])
+    end
+  end
+
+  return H
+end
+
+"""
+    instantiate_MPO(H_template::MPO, sym::SymbolicMPO, coefficients; checkflux=false) -> MPO
+
+Copy a compatible MPO template and evaluate `sym` into the copied layout.
+
+The template fixes the tensor block layout, so this overload does not accept
+`splitblocks`.
+"""
+function instantiate_MPO(
+  H_template::MPO, sym::SymbolicMPO, coefficients::AbstractVector; checkflux::Bool=false
+)::MPO
+  return instantiate_MPO!(deepcopy(H_template), sym, coefficients; checkflux)
+end
+
 function _check_MPO_symbolic_kwargs(kwargs)::Nothing
   if haskey(kwargs, :alg)
     throw(
diff --git a/test/symbolic-local-matrix-tests.jl b/test/symbolic-local-matrix-tests.jl
index 6117eee..72cdf65 100644
--- a/test/symbolic-local-matrix-tests.jl
+++ b/test/symbolic-local-matrix-tests.jl
@@ -394,7 +394,7 @@ function mpo_relative_error(A::MPO, B::MPO)::Float64
   AmB = add(A, -B; alg="directsum")
   numerator = real(inner(AmB, AmB))
   denominator = real(inner(A, A))
-  return iszero(denominator) ? numerator : numerator / denominator
+  return numerator / denominator
 end
 
 function test_mpos_approx_equal(A::MPO, B::MPO; tol::Real=1e-10)::Nothing
@@ -518,6 +518,57 @@ function test_symbolic_mpo_qn_fresh_instantiation()::Nothing
   return nothing
 end
 
+function test_symbolic_mpo_inplace_instantiation()::Nothing
+  os, sites = qn_number_symbolic_fixture()
+  sym = MPO_symbolic(os, sites)
+
+  coefficients1 = [0.8, -1.1, 0.3]
+  coefficients2 = [-0.2, 1.4, 0.6]
+
+  H = instantiate_MPO(sym, coefficients1; splitblocks=true)
+  original_objectid = objectid(H)
+  original_linkdims = linkdims(H)
+
+  result = instantiate_MPO!(H, sym, coefficients2; checkflux=true)
+  expected = instantiate_MPO(sym, coefficients2; splitblocks=true, checkflux=true)
+
+  @test result === H
+  @test objectid(H) == original_objectid
+  @test linkdims(H) == original_linkdims
+  test_mpos_approx_equal(H, expected)
+
+  @test_throws MethodError instantiate_MPO!(H, sym, coefficients2; splitblocks=false)
+
+  return nothing
+end
+
+function test_symbolic_mpo_template_instantiation()::Nothing
+  os, sites = qn_number_symbolic_fixture()
+  sym = MPO_symbolic(os, sites)
+
+  coefficients1 = [0.8, -1.1, 0.3]
+  coefficients2 = [-0.2, 1.4, 0.6]
+
+  H_template = instantiate_MPO(sym, coefficients1; splitblocks=false)
+  H_template_expected = instantiate_MPO(sym, coefficients1; splitblocks=false)
+  H = instantiate_MPO(H_template, sym, coefficients2; checkflux=true)
+  expected = instantiate_MPO(sym, coefficients2; splitblocks=false, checkflux=true)
+
+  @test H !== H_template
+  @test linkdims(H) == linkdims(H_template)
+  test_mpos_approx_equal(H, expected)
+  test_mpos_approx_equal(H_template, H_template_expected)
+
+  @test_throws MethodError instantiate_MPO(H_template, sym, coefficients2; splitblocks=true)
+
+  bad_os, bad_sites = qn_number_symbolic_fixture()
+  bad_sym = MPO_symbolic(bad_os, bad_sites)
+  H_bad = instantiate_MPO(bad_sym, coefficients1; splitblocks=false)
+  @test_throws Exception instantiate_MPO(H_bad, sym, coefficients2)
+
+  return nothing
+end
+
 function test_symbolic_mpo_uses_negative_local_op_id_for_jw_terms()::Nothing
   sites = siteinds("Fermion", 2)
   op_cache_vec = to_OpCacheVec(sites, [["I", "C", "Cdag", "N"] for _ in eachindex(sites)])
@@ -562,5 +613,9 @@ end
 
   test_symbolic_mpo_qn_fresh_instantiation()
 
+  test_symbolic_mpo_inplace_instantiation()
+
+  test_symbolic_mpo_template_instantiation()
+
   test_symbolic_mpo_uses_negative_local_op_id_for_jw_terms()
 end

From 25b6ee24c7be5c32dd4b5fd1fc36f2c9ddc405fc Mon Sep 17 00:00:00 2001
From: Ben Corbett <corbett5@llnl.gov>
Date: Mon, 27 Apr 2026 14:56:55 -0600
Subject: [PATCH 09/16] Step 8

---
 src/MPOConstruction.jl              |  32 ++-
 test/runtests.jl                    |   1 +
 test/symbolic-local-matrix-tests.jl | 308 +++------------------------
 test/symbolic-mpo-tests.jl          | 310 ++++++++++++++++++++++++++++
 4 files changed, 365 insertions(+), 286 deletions(-)
 create mode 100644 test/symbolic-mpo-tests.jl

diff --git a/src/MPOConstruction.jl b/src/MPOConstruction.jl
index 02e6911..6a17302 100644
--- a/src/MPOConstruction.jl
+++ b/src/MPOConstruction.jl
@@ -314,6 +314,31 @@ function _fill_symbolic_mpo_tensor_with_boundary_links!(
   return tensor
 end
 
+function _check_symbolic_mpo_template(H::MPO, sym::SymbolicMPO)::Nothing
+  if length(H) != length(sym.sites)
+    throw(
+      ArgumentError(
+        "Template MPO length $(length(H)) is incompatible with SymbolicMPO length $(length(sym.sites)).",
+      ),
+    )
+  end
+
+  for n in eachindex(sym.sites)
+    if !hasinds(H[n], prime(sym.sites[n]), dag(sym.sites[n]))
+      throw(ArgumentError("Template MPO site index at site $n is incompatible with SymbolicMPO."))
+    end
+  end
+
+  for n in 2:(length(sym.llinks) - 1)
+    template_link = linkind(H, n - 1)
+    if isnothing(template_link) || template_link != sym.llinks[n]
+      throw(ArgumentError("Template MPO link index at bond $(n - 1) is incompatible with SymbolicMPO."))
+    end
+  end
+
+  return nothing
+end
+
 """
     instantiate_MPO!(H::MPO, sym::SymbolicMPO, coefficients; checkflux=false) -> MPO
 
@@ -330,12 +355,7 @@ function instantiate_MPO!(
   H::MPO, sym::SymbolicMPO, coefficients::AbstractVector; checkflux::Bool=false
 )::MPO
   _check_symbolic_coefficients(sym, coefficients)
-
-  if length(H) != length(sym.sites)
-    throw(ArgumentError(
-      "Template MPO length $(length(H)) is incompatible with SymbolicMPO length $(length(sym.sites)).",
-    ))
-  end
+  _check_symbolic_mpo_template(H, sym)
 
   _add_symbolic_mpo_boundary_links!(H, sym)
 
diff --git a/test/runtests.jl b/test/runtests.jl
index fd80071..61c9eb7 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,6 +1,7 @@
 include("tensor-construction-tests.jl")
 include("ops-tests.jl")
 include("symbolic-local-matrix-tests.jl")
+include("symbolic-mpo-tests.jl")
 include("large-graph-tests.jl")
 include("test-MPOConstruction.jl")
 include("mponew_h1_accuracy.jl")
diff --git a/test/symbolic-local-matrix-tests.jl b/test/symbolic-local-matrix-tests.jl
index 72cdf65..21b66bc 100644
--- a/test/symbolic-local-matrix-tests.jl
+++ b/test/symbolic-local-matrix-tests.jl
@@ -1,8 +1,6 @@
 using ITensorMPOConstruction
 using ITensorMPOConstruction:
   MPOGraph,
-  MPO_symbolic,
-  SymbolicMPO,
   SymbolicLocalMatrix,
   _append_symbolic_local_matrix_term!,
   _evaluate_symbolic_local_matrix,
@@ -121,6 +119,19 @@ function test_internalize_symbolic_ids()::Nothing
   return nothing
 end
 
+function test_negative_internal_symbolic_ids_remain_supported()::Nothing
+  coefficients = [7.0]
+
+  @test _weight_value(-2, coefficients) == -coefficients[1]
+
+  terms = SymbolicLocalMatrix{Int}()
+  _append_symbolic_local_matrix_term!(terms, -2, 2)
+  @test terms == [(-2, 2)]
+  @test _max_user_label(terms) == 1
+
+  return nothing
+end
+
 function symbolic_basis_rewrite_cache_vecs(factor)
   identity_matrix = [1.0 0.0; 0.0 1.0]
   x_matrix = [0.0 1.0; 1.0 0.0]
@@ -334,288 +345,25 @@ function test_symbolic_mpo_graph_cancels_opposite_signed_ids()::Nothing
   return nothing
 end
 
-function transverse_field_ising_symbolic_fixture()
-  sites = siteinds("Qubit", 2)
-  op_cache_vec = to_OpCacheVec(sites, [["I", "X", "Z"] for _ in eachindex(sites)])
-
-  X(n) = OpID(2, n)
-  Z(n) = OpID(3, n)
-
-  os = OpIDSum{2,Int,Int}(3, op_cache_vec)
-  add!(os, 1, X(1), X(2))
-  add!(os, 2, Z(1))
-  add!(os, 3, Z(2))
-
-  return os, sites
-end
-
-function transverse_field_ising_numeric_opidsum(sites, coefficients)
-  op_cache_vec = to_OpCacheVec(sites, [["I", "X", "Z"] for _ in eachindex(sites)])
-
-  X(n) = OpID(2, n)
-  Z(n) = OpID(3, n)
-
-  os = OpIDSum{2,eltype(coefficients),Int}(3, op_cache_vec)
-  add!(os, coefficients[1], X(1), X(2))
-  add!(os, coefficients[2], Z(1))
-  add!(os, coefficients[3], Z(2))
-
-  return os
-end
-
-function qn_number_symbolic_fixture()
-  sites = siteinds("Fermion", 2; conserve_qns=true)
-  op_cache_vec = to_OpCacheVec(sites, [["I", "N"] for _ in eachindex(sites)])
-
-  Nop(n) = OpID(2, n)
-
-  os = OpIDSum{2,Int,Int}(3, op_cache_vec)
-  add!(os, 1, Nop(1))
-  add!(os, 2, Nop(2))
-  add!(os, 3, Nop(1), Nop(2))
-
-  return os, sites
-end
-
-function qn_number_numeric_opidsum(sites, coefficients)
-  op_cache_vec = to_OpCacheVec(sites, [["I", "N"] for _ in eachindex(sites)])
-
-  Nop(n) = OpID(2, n)
-
-  os = OpIDSum{2,eltype(coefficients),Int}(3, op_cache_vec)
-  add!(os, coefficients[1], Nop(1))
-  add!(os, coefficients[2], Nop(2))
-  add!(os, coefficients[3], Nop(1), Nop(2))
-
-  return os
-end
-
-function mpo_relative_error(A::MPO, B::MPO)::Float64
-  AmB = add(A, -B; alg="directsum")
-  numerator = real(inner(AmB, AmB))
-  denominator = real(inner(A, A))
-  return numerator / denominator
-end
-
-function test_mpos_approx_equal(A::MPO, B::MPO; tol::Real=1e-10)::Nothing
-  relative_error = mpo_relative_error(A, B)
-  @test relative_error < tol
-  return nothing
-end
-
-function symbolic_mpo_terms(sym::SymbolicMPO)::Vector{Tuple{Int,Int}}
-  terms = Tuple{Int,Int}[]
-  for site_matrices in sym.block_sparse_matrices
-    for matrix_of_component in site_matrices
-      for block in matrix_of_component
-        for block_terms in values(block)
-          append!(terms, block_terms)
-        end
-      end
-    end
+@testset "SymbolicLocalMatrix" begin
+  @testset "symbolic local matrix terms" begin
+    test_symbolic_local_matrix_operations()
   end
-  return terms
-end
-
-function test_symbolic_mpo_construction_metadata()::Nothing
-  os, sites = transverse_field_ising_symbolic_fixture()
 
-  sym = MPO_symbolic(os, sites)
-
-  @test sym isa SymbolicMPO
-  @test length(sym.offsets) == length(sites)
-  @test length(sym.block_sparse_matrices) == length(sites)
-  @test length(sym.llinks) == length(sites) + 1
-  @test sym.max_user_label == 3
-  @test sym.op_cache_vec === os.op_cache_vec
-  @test !isempty(symbolic_mpo_terms(sym))
-
-  return nothing
-end
-
-function test_symbolic_mpo_rejects_unsupported_public_inputs()::Nothing
-  os, sites = transverse_field_ising_symbolic_fixture()
-  @test_throws ArgumentError MPO_symbolic(os, sites; alg="QR")
-
-  noninteger_os = OpIDSum{2,Float64,Int}(1, os.op_cache_vec)
-  add!(noninteger_os, 1.0, OpID(2, 1), OpID(2, 2))
-  @test_throws ArgumentError MPO_symbolic(noninteger_os, sites)
-
-  op_sum = OpSum()
-  op_sum += "X", 1
-  @test_throws ArgumentError MPO_symbolic(op_sum, sites)
-
-  return nothing
-end
-
-function test_symbolic_mpo_fresh_instantiation_matches_numeric()::Nothing
-  os, sites = transverse_field_ising_symbolic_fixture()
-  sym = MPO_symbolic(os, sites)
-
-  for coefficients in ([1.25, -0.5, 0.75], [-2.0, 0.25, 1.5])
-    symbolic_mpo = instantiate_MPO(sym, coefficients; splitblocks=true)
-    numeric_mpo = MPO_new(
-      transverse_field_ising_numeric_opidsum(sites, coefficients),
-      sites;
-      alg="VC",
-      splitblocks=true,
-      checkflux=false,
-    )
-
-    test_mpos_approx_equal(symbolic_mpo, numeric_mpo)
+  @testset "internal symbolic ids" begin
+    test_internalize_symbolic_ids()
+    test_negative_internal_symbolic_ids_remain_supported()
   end
 
-  return nothing
-end
-
-function test_symbolic_mpo_duplicate_labels_instantiate_like_numeric()::Nothing
-  sites, op_cache_vec = two_site_qubit_symbolic_fixture()
-
-  X(n) = OpID(2, n)
-
-  symbolic_os = OpIDSum{2,Int,Int}(2, op_cache_vec)
-  add!(symbolic_os, 1, X(1), X(2))
-  add!(symbolic_os, 2, X(1), X(2))
-
-  sym = MPO_symbolic(symbolic_os, sites)
-  coefficients = [2.0, 3.0]
-  symbolic_mpo = instantiate_MPO(sym, coefficients; splitblocks=true)
-
-  numeric_os = OpIDSum{2,Float64,Int}(2, op_cache_vec)
-  add!(numeric_os, coefficients[1], X(1), X(2))
-  add!(numeric_os, coefficients[2], X(1), X(2))
-  numeric_mpo = MPO_new(numeric_os, sites; alg="VC", splitblocks=true, checkflux=false)
-
-  test_mpos_approx_equal(symbolic_mpo, numeric_mpo)
-
-  return nothing
-end
-
-function test_symbolic_mpo_qn_fresh_instantiation()::Nothing
-  os, sites = qn_number_symbolic_fixture()
-  sym = MPO_symbolic(os, sites)
-
-  coefficients = [0.8, -1.1, 0.3]
-  for splitblocks in (false, true)
-    symbolic_mpo = instantiate_MPO(sym, coefficients; splitblocks, checkflux=true)
-    numeric_mpo = MPO_new(
-      qn_number_numeric_opidsum(sites, coefficients),
-      sites;
-      alg="VC",
-      splitblocks,
-      checkflux=true,
-    )
-
-    test_mpos_approx_equal(symbolic_mpo, numeric_mpo)
+  @testset "basis rewrite" begin
+    test_symbolic_basis_rewrite_sign_factors()
+    test_symbolic_basis_rewrite_rejects_unsupported_factors()
   end
 
-  nonzero_mpo = instantiate_MPO(sym, [1.0, 2.0, 3.0]; splitblocks=true)
-  zero_mpo = instantiate_MPO(sym, [0.0, 0.0, 0.0]; splitblocks=true)
-
-  @test linkdims(zero_mpo) == linkdims(nonzero_mpo)
-  @test_throws ArgumentError instantiate_MPO(sym, [1.0, 2.0])
-
-  return nothing
-end
-
-function test_symbolic_mpo_inplace_instantiation()::Nothing
-  os, sites = qn_number_symbolic_fixture()
-  sym = MPO_symbolic(os, sites)
-
-  coefficients1 = [0.8, -1.1, 0.3]
-  coefficients2 = [-0.2, 1.4, 0.6]
-
-  H = instantiate_MPO(sym, coefficients1; splitblocks=true)
-  original_objectid = objectid(H)
-  original_linkdims = linkdims(H)
-
-  result = instantiate_MPO!(H, sym, coefficients2; checkflux=true)
-  expected = instantiate_MPO(sym, coefficients2; splitblocks=true, checkflux=true)
-
-  @test result === H
-  @test objectid(H) == original_objectid
-  @test linkdims(H) == original_linkdims
-  test_mpos_approx_equal(H, expected)
-
-  @test_throws MethodError instantiate_MPO!(H, sym, coefficients2; splitblocks=false)
-
-  return nothing
-end
-
-function test_symbolic_mpo_template_instantiation()::Nothing
-  os, sites = qn_number_symbolic_fixture()
-  sym = MPO_symbolic(os, sites)
-
-  coefficients1 = [0.8, -1.1, 0.3]
-  coefficients2 = [-0.2, 1.4, 0.6]
-
-  H_template = instantiate_MPO(sym, coefficients1; splitblocks=false)
-  H_template_expected = instantiate_MPO(sym, coefficients1; splitblocks=false)
-  H = instantiate_MPO(H_template, sym, coefficients2; checkflux=true)
-  expected = instantiate_MPO(sym, coefficients2; splitblocks=false, checkflux=true)
-
-  @test H !== H_template
-  @test linkdims(H) == linkdims(H_template)
-  test_mpos_approx_equal(H, expected)
-  test_mpos_approx_equal(H_template, H_template_expected)
-
-  @test_throws MethodError instantiate_MPO(H_template, sym, coefficients2; splitblocks=true)
-
-  bad_os, bad_sites = qn_number_symbolic_fixture()
-  bad_sym = MPO_symbolic(bad_os, bad_sites)
-  H_bad = instantiate_MPO(bad_sym, coefficients1; splitblocks=false)
-  @test_throws Exception instantiate_MPO(H_bad, sym, coefficients2)
-
-  return nothing
-end
-
-function test_symbolic_mpo_uses_negative_local_op_id_for_jw_terms()::Nothing
-  sites = siteinds("Fermion", 2)
-  op_cache_vec = to_OpCacheVec(sites, [["I", "C", "Cdag", "N"] for _ in eachindex(sites)])
-
-  C(n) = OpID(2, n)
-  Cdag(n) = OpID(3, n)
-
-  os = OpIDSum{2,Int,Int}(1, op_cache_vec)
-  add!(os, 1, Cdag(1), C(2))
-
-  sym = MPO_symbolic(os, sites)
-  terms = symbolic_mpo_terms(sym)
-  @test any(term -> term[2] < 0, terms)
-
-  return nothing
-end
-
-@testset "SymbolicLocalMatrix" begin
-  test_symbolic_local_matrix_operations()
-
-  test_internalize_symbolic_ids()
-
-  test_symbolic_basis_rewrite_sign_factors()
-
-  test_symbolic_basis_rewrite_rejects_unsupported_factors()
-
-  test_numeric_mpo_graph_still_compacts_duplicate_terms()
-
-  test_symbolic_mpo_graph_preserves_duplicate_signed_ids()
-
-  test_symbolic_mpo_graph_preserves_duplicate_label_multiplicity()
-
-  test_symbolic_mpo_graph_cancels_opposite_signed_ids()
-
-  test_symbolic_mpo_construction_metadata()
-
-  test_symbolic_mpo_rejects_unsupported_public_inputs()
-
-  test_symbolic_mpo_fresh_instantiation_matches_numeric()
-
-  test_symbolic_mpo_duplicate_labels_instantiate_like_numeric()
-
-  test_symbolic_mpo_qn_fresh_instantiation()
-
-  test_symbolic_mpo_inplace_instantiation()
-
-  test_symbolic_mpo_template_instantiation()
-
-  test_symbolic_mpo_uses_negative_local_op_id_for_jw_terms()
+  @testset "symbolic graph duplicate terms" begin
+    test_numeric_mpo_graph_still_compacts_duplicate_terms()
+    test_symbolic_mpo_graph_preserves_duplicate_signed_ids()
+    test_symbolic_mpo_graph_preserves_duplicate_label_multiplicity()
+    test_symbolic_mpo_graph_cancels_opposite_signed_ids()
+  end
 end
diff --git a/test/symbolic-mpo-tests.jl b/test/symbolic-mpo-tests.jl
new file mode 100644
index 0000000..a4c5625
--- /dev/null
+++ b/test/symbolic-mpo-tests.jl
@@ -0,0 +1,310 @@
+using ITensorMPOConstruction
+using ITensorMPOConstruction: MPO_symbolic, SymbolicMPO
+using ITensors
+using ITensorMPS
+using Test
+
+function two_site_qubit_symbolic_mpo_fixture()
+  sites = siteinds("Qubit", 2)
+  op_cache_vec = to_OpCacheVec(sites, [["I", "X"] for _ in eachindex(sites)])
+  return sites, op_cache_vec
+end
+
+function transverse_field_ising_symbolic_fixture()
+  sites = siteinds("Qubit", 2)
+  op_cache_vec = to_OpCacheVec(sites, [["I", "X", "Z"] for _ in eachindex(sites)])
+
+  X(n) = OpID(2, n)
+  Z(n) = OpID(3, n)
+
+  os = OpIDSum{2,Int,Int}(3, op_cache_vec)
+  add!(os, 1, X(1), X(2))
+  add!(os, 2, Z(1))
+  add!(os, 3, Z(2))
+
+  return os, sites
+end
+
+function transverse_field_ising_numeric_opidsum(sites, coefficients)
+  op_cache_vec = to_OpCacheVec(sites, [["I", "X", "Z"] for _ in eachindex(sites)])
+
+  X(n) = OpID(2, n)
+  Z(n) = OpID(3, n)
+
+  os = OpIDSum{2,eltype(coefficients),Int}(3, op_cache_vec)
+  add!(os, coefficients[1], X(1), X(2))
+  add!(os, coefficients[2], Z(1))
+  add!(os, coefficients[3], Z(2))
+
+  return os
+end
+
+function qn_number_symbolic_fixture()
+  sites = siteinds("Fermion", 2; conserve_qns=true)
+  op_cache_vec = to_OpCacheVec(sites, [["I", "N"] for _ in eachindex(sites)])
+
+  Nop(n) = OpID(2, n)
+
+  os = OpIDSum{2,Int,Int}(3, op_cache_vec)
+  add!(os, 1, Nop(1))
+  add!(os, 2, Nop(2))
+  add!(os, 3, Nop(1), Nop(2))
+
+  return os, sites
+end
+
+function qn_number_numeric_opidsum(sites, coefficients)
+  op_cache_vec = to_OpCacheVec(sites, [["I", "N"] for _ in eachindex(sites)])
+
+  Nop(n) = OpID(2, n)
+
+  os = OpIDSum{2,eltype(coefficients),Int}(3, op_cache_vec)
+  add!(os, coefficients[1], Nop(1))
+  add!(os, coefficients[2], Nop(2))
+  add!(os, coefficients[3], Nop(1), Nop(2))
+
+  return os
+end
+
+function mpo_relative_error(A::MPO, B::MPO)::Float64
+  AmB = add(A, -B; alg="directsum")
+  numerator = real(inner(AmB, AmB))
+  denominator = real(inner(A, A))
+  return numerator / denominator
+end
+
+function test_mpos_approx_equal(A::MPO, B::MPO; tol::Real=1e-10)::Nothing
+  relative_error = mpo_relative_error(A, B)
+  @test relative_error < tol
+  return nothing
+end
+
+function symbolic_mpo_terms(sym::SymbolicMPO)::Vector{Tuple{Int,Int}}
+  terms = Tuple{Int,Int}[]
+  for site_matrices in sym.block_sparse_matrices
+    for matrix_of_component in site_matrices
+      for block in matrix_of_component
+        for block_terms in values(block)
+          append!(terms, block_terms)
+        end
+      end
+    end
+  end
+  return terms
+end
+
+function test_symbolic_mpo_construction_metadata()::Nothing
+  os, sites = transverse_field_ising_symbolic_fixture()
+
+  sym = MPO_symbolic(os, sites)
+
+  @test sym isa SymbolicMPO
+  @test length(sym.offsets) == length(sites)
+  @test length(sym.block_sparse_matrices) == length(sites)
+  @test length(sym.llinks) == length(sites) + 1
+  @test sym.max_user_label == 3
+  @test sym.op_cache_vec === os.op_cache_vec
+  @test !isempty(symbolic_mpo_terms(sym))
+
+  return nothing
+end
+
+function test_symbolic_mpo_rejects_unsupported_public_inputs()::Nothing
+  os, sites = transverse_field_ising_symbolic_fixture()
+  @test_throws ArgumentError MPO_symbolic(os, sites; alg="QR")
+
+  noninteger_os = OpIDSum{2,Float64,Int}(1, os.op_cache_vec)
+  add!(noninteger_os, 1.0, OpID(2, 1), OpID(2, 2))
+  @test_throws ArgumentError MPO_symbolic(noninteger_os, sites)
+
+  op_sum = OpSum()
+  op_sum += "X", 1
+  @test_throws ArgumentError MPO_symbolic(op_sum, sites)
+
+  return nothing
+end
+
+function test_symbolic_mpo_fresh_instantiation_matches_numeric()::Nothing
+  os, sites = transverse_field_ising_symbolic_fixture()
+  sym = MPO_symbolic(os, sites)
+
+  for coefficients in ([1.25, -0.5, 0.75], [-2.0, 0.25, 1.5])
+    symbolic_mpo = instantiate_MPO(sym, coefficients; splitblocks=true)
+    numeric_mpo = MPO_new(
+      transverse_field_ising_numeric_opidsum(sites, coefficients),
+      sites;
+      alg="VC",
+      splitblocks=true,
+      checkflux=false,
+    )
+
+    test_mpos_approx_equal(symbolic_mpo, numeric_mpo)
+  end
+
+  return nothing
+end
+
+function test_symbolic_mpo_duplicate_labels_instantiate_like_numeric()::Nothing
+  sites, op_cache_vec = two_site_qubit_symbolic_mpo_fixture()
+
+  X(n) = OpID(2, n)
+
+  symbolic_os = OpIDSum{2,Int,Int}(2, op_cache_vec)
+  add!(symbolic_os, 1, X(1), X(2))
+  add!(symbolic_os, 2, X(1), X(2))
+
+  sym = MPO_symbolic(symbolic_os, sites)
+  coefficients = [2.0, 3.0]
+  symbolic_mpo = instantiate_MPO(sym, coefficients; splitblocks=true)
+
+  numeric_os = OpIDSum{2,Float64,Int}(2, op_cache_vec)
+  add!(numeric_os, coefficients[1], X(1), X(2))
+  add!(numeric_os, coefficients[2], X(1), X(2))
+  numeric_mpo = MPO_new(numeric_os, sites; alg="VC", splitblocks=true, checkflux=false)
+
+  test_mpos_approx_equal(symbolic_mpo, numeric_mpo)
+
+  return nothing
+end
+
+function test_symbolic_mpo_qn_fresh_instantiation()::Nothing
+  os, sites = qn_number_symbolic_fixture()
+  sym = MPO_symbolic(os, sites)
+
+  coefficients = [0.8, -1.1, 0.3]
+  for splitblocks in (false, true)
+    symbolic_mpo = instantiate_MPO(sym, coefficients; splitblocks, checkflux=true)
+    numeric_mpo = MPO_new(
+      qn_number_numeric_opidsum(sites, coefficients),
+      sites;
+      alg="VC",
+      splitblocks,
+      checkflux=true,
+    )
+
+    test_mpos_approx_equal(symbolic_mpo, numeric_mpo)
+  end
+
+  nonzero_mpo = instantiate_MPO(sym, [1.0, 2.0, 3.0]; splitblocks=true)
+  zero_mpo = instantiate_MPO(sym, [0.0, 0.0, 0.0]; splitblocks=true)
+
+  @test linkdims(zero_mpo) == linkdims(nonzero_mpo)
+  @test_throws ArgumentError instantiate_MPO(sym, [1.0, 2.0])
+
+  return nothing
+end
+
+function test_symbolic_mpo_missing_coefficients_rejected()::Nothing
+  os, sites = transverse_field_ising_symbolic_fixture()
+  sym = MPO_symbolic(os, sites)
+  coefficients = [1.0, 2.0, 3.0]
+  short_coefficients = [1.0, 2.0]
+
+  H = instantiate_MPO(sym, coefficients; splitblocks=true)
+
+  @test_throws ArgumentError instantiate_MPO(sym, short_coefficients; splitblocks=true)
+  @test_throws ArgumentError instantiate_MPO!(H, sym, short_coefficients)
+  @test_throws ArgumentError instantiate_MPO(H, sym, short_coefficients)
+
+  return nothing
+end
+
+function test_symbolic_mpo_inplace_instantiation()::Nothing
+  os, sites = qn_number_symbolic_fixture()
+  sym = MPO_symbolic(os, sites)
+
+  coefficients1 = [0.8, -1.1, 0.3]
+  coefficients2 = [-0.2, 1.4, 0.6]
+
+  H = instantiate_MPO(sym, coefficients1; splitblocks=true)
+  original_objectid = objectid(H)
+  original_linkdims = linkdims(H)
+
+  result = instantiate_MPO!(H, sym, coefficients2; checkflux=true)
+  expected = instantiate_MPO(sym, coefficients2; splitblocks=true, checkflux=true)
+
+  @test result === H
+  @test objectid(H) == original_objectid
+  @test linkdims(H) == original_linkdims
+  test_mpos_approx_equal(H, expected)
+
+  @test_throws MethodError instantiate_MPO!(H, sym, coefficients2; splitblocks=false)
+
+  return nothing
+end
+
+function test_symbolic_mpo_template_instantiation()::Nothing
+  os, sites = qn_number_symbolic_fixture()
+  sym = MPO_symbolic(os, sites)
+
+  coefficients1 = [0.8, -1.1, 0.3]
+  coefficients2 = [-0.2, 1.4, 0.6]
+
+  H_template = instantiate_MPO(sym, coefficients1; splitblocks=false)
+  H_template_expected = instantiate_MPO(sym, coefficients1; splitblocks=false)
+  H = instantiate_MPO(H_template, sym, coefficients2; checkflux=true)
+  expected = instantiate_MPO(sym, coefficients2; splitblocks=false, checkflux=true)
+
+  @test H !== H_template
+  @test linkdims(H) == linkdims(H_template)
+  test_mpos_approx_equal(H, expected)
+  test_mpos_approx_equal(H_template, H_template_expected)
+
+  @test_throws MethodError instantiate_MPO(H_template, sym, coefficients2; splitblocks=true)
+
+  return nothing
+end
+
+function test_symbolic_mpo_incompatible_template_rejected()::Nothing
+  os, sites = qn_number_symbolic_fixture()
+  sym = MPO_symbolic(os, sites)
+  coefficients = [0.8, -1.1, 0.3]
+
+  bad_os, bad_sites = qn_number_symbolic_fixture()
+  bad_sym = MPO_symbolic(bad_os, bad_sites)
+  H_bad = instantiate_MPO(bad_sym, coefficients; splitblocks=false)
+  @test_throws ArgumentError instantiate_MPO(H_bad, sym, coefficients)
+
+  return nothing
+end
+
+function test_symbolic_mpo_uses_negative_local_op_id_for_jw_terms()::Nothing
+  sites = siteinds("Fermion", 2)
+  op_cache_vec = to_OpCacheVec(sites, [["I", "C", "Cdag", "N"] for _ in eachindex(sites)])
+
+  C(n) = OpID(2, n)
+  Cdag(n) = OpID(3, n)
+
+  os = OpIDSum{2,Int,Int}(1, op_cache_vec)
+  add!(os, 1, Cdag(1), C(2))
+
+  sym = MPO_symbolic(os, sites)
+  terms = symbolic_mpo_terms(sym)
+  @test any(term -> term[2] < 0, terms)
+
+  return nothing
+end
+
+@testset "SymbolicMPO" begin
+  @testset "fresh instantiation" begin
+    test_symbolic_mpo_construction_metadata()
+    test_symbolic_mpo_fresh_instantiation_matches_numeric()
+    test_symbolic_mpo_duplicate_labels_instantiate_like_numeric()
+    test_symbolic_mpo_qn_fresh_instantiation()
+    test_symbolic_mpo_uses_negative_local_op_id_for_jw_terms()
+  end
+
+  @testset "public API rejections" begin
+    test_symbolic_mpo_rejects_unsupported_public_inputs()
+    test_symbolic_mpo_missing_coefficients_rejected()
+    test_symbolic_mpo_incompatible_template_rejected()
+  end
+
+  @testset "in-place instantiation" begin
+    test_symbolic_mpo_inplace_instantiation()
+  end
+
+  @testset "template-assisted instantiation" begin
+    test_symbolic_mpo_template_instantiation()
+  end
+end

From 27f9824785dc185cb12373ec8c441429b09da6a1 Mon Sep 17 00:00:00 2001
From: Ben Corbett <corbett5@llnl.gov>
Date: Mon, 27 Apr 2026 15:55:18 -0600
Subject: [PATCH 10/16] added TFIM example.

---
 docs/make.jl                                |   2 +
 docs/src/documentation/MPO_new.md           |   9 ++
 docs/src/documentation/unexported.md        |   3 +
 examples/symbolic-transverse-field-ising.jl | 108 ++++++++++++++++++++
 src/MPOConstruction.jl                      |   2 +
 5 files changed, 124 insertions(+)
 create mode 100644 examples/symbolic-transverse-field-ising.jl

diff --git a/docs/make.jl b/docs/make.jl
index e82688b..bc591bf 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -17,6 +17,7 @@ examples = [
   "./examples/fermi-hubbard-ks.jl",
   "./examples/haldane-shastry.jl",
   "./examples/electronic-structure.jl",
+  "./examples/symbolic-transverse-field-ising.jl",
   "./examples/fermi-hubbard-tc.jl",
 ]
 
@@ -53,6 +54,7 @@ makedocs(;
       "Real Space Fermi-Hubbard" => "examples/fermi-hubbard-rs.md",
       "Momentum Space Fermi-Hubbard" => "examples/fermi-hubbard-ks.md",
       "Electronic Structure" => "examples/electronic-structure.md",
+      "Symbolic Transverse Field Ising" => "examples/symbolic-transverse-field-ising.md",
       "Haldane-Shastry and Truncation" => "examples/haldane-shastry.md",
       "Challenge Problem" => "examples/fermi-hubbard-tc.md",
     ],
diff --git a/docs/src/documentation/MPO_new.md b/docs/src/documentation/MPO_new.md
index 49ddafc..509c3a6 100644
--- a/docs/src/documentation/MPO_new.md
+++ b/docs/src/documentation/MPO_new.md
@@ -9,11 +9,20 @@ The main exported function is `MPO_new`, which takes an `OpSum` and transforms i
 
 ```@docs
 MPO_new
+SymbolicMPO
+MPO_symbolic
 resume_MPO_construction!
 instantiate_MPO
+instantiate_MPO!
 sparsity
 block2_nnz
 ```
 
+## Symbolic MPO construction
+`MPO_symbolic` constructs the vertex-cover MPO structure once from an integer-labeled
+`OpIDSum`. The labels are positive one-based coefficient ids, and
+`instantiate_MPO` or `instantiate_MPO!` substitutes numeric values later without
+rebuilding the graph.
+
 ## A note on truncation
 `MPO_new` is designed to construct _numerically exact_ MPOs; the tolerance parameter should only be used to adjust the definition of "numerically exact", not to perform truncation. If truncation is desired, truncate the resulting MPO with `ITensorMPS.truncate!`. See [Haldane-Shastry and truncation](../examples/haldane-shastry.md) for an example.
diff --git a/docs/src/documentation/unexported.md b/docs/src/documentation/unexported.md
index 9f045f3..7c9c2a4 100644
--- a/docs/src/documentation/unexported.md
+++ b/docs/src/documentation/unexported.md
@@ -13,6 +13,9 @@ ITensorMPOConstruction.rewrite_in_operator_basis!
 ITensorMPOConstruction.op_sum_to_opID_sum
 ITensorMPOConstruction.check_os_for_errors
 ITensorMPOConstruction.prepare_opID_sum!
+ITensorMPOConstruction.SymbolicLocalMatrix
+ITensorMPOConstruction.SymbolicBlockSparseMatrix
+ITensorMPOConstruction.internalize_symbolic_ids!
 ITensorMPOConstruction.get_onsite_op
 ```
 
diff --git a/examples/symbolic-transverse-field-ising.jl b/examples/symbolic-transverse-field-ising.jl
new file mode 100644
index 0000000..32110e3
--- /dev/null
+++ b/examples/symbolic-transverse-field-ising.jl
@@ -0,0 +1,108 @@
+# # Symbolic all-to-all transverse-field Ising model
+#
+# The all-to-all transverse-field Ising model is
+# ```math
+# H = \sum_{i < j} J_{ij} Z_i Z_j + h \sum_i X_i .
+# ```
+# If the connectivity changes many times while the operator structure stays
+# fixed, the symbolic MPO path can construct the MPO once and instantiate
+# it repeatedly for different $J_{ij}$ and $h$.
+
+using ITensors, ITensorMPS, ITensorMPOConstruction
+using Random
+
+function num_tfim_couplings(N::Int)::Int
+  return (N * (N - 1)) ÷ 2
+end
+
+function all_to_all_tfim_symbolic(N::Int)::SymbolicMPO
+  sites = siteinds("Qubit", N)
+  op_cache_vec = to_OpCacheVec(sites, [["I", "X", "Z"] for _ in 1:N])
+
+  X(n::Int) = OpID(2, n)
+  Z(n::Int) = OpID(3, n)
+
+  num_couplings = num_tfim_couplings(N)
+  field_label = num_couplings + 1
+  os = OpIDSum{2,Int,Int}(num_couplings + N, op_cache_vec)
+
+  label = 1
+  for i in 1:N
+    for j in (i + 1):N
+      add!(os, label, Z(i), Z(j))
+      label += 1
+    end
+  end
+
+  for i in 1:N
+    add!(os, field_label, X(i))
+  end
+
+  return MPO_symbolic(os, sites)
+end
+
+function all_to_all_tfim_numeric(
+  sites::Vector{<:Index}, coupling_weights::AbstractVector, h::Real
+)::MPO
+  N = length(sites)
+  num_couplings = num_tfim_couplings(N)
+
+  op_cache_vec = to_OpCacheVec(sites, [["I", "X", "Z"] for _ in 1:N])
+
+  X(n::Int) = OpID(2, n)
+  Z(n::Int) = OpID(3, n)
+
+  os = OpIDSum{2,Float64,Int}(num_couplings + N, op_cache_vec)
+
+  label = 1
+  for i in 1:N
+    for j in (i + 1):N
+      add!(os, coupling_weights[label], Z(i), Z(j))
+      label += 1
+    end
+  end
+
+  for i in 1:N
+    add!(os, h, X(i))
+  end
+
+  return MPO_new(os, sites; alg="VC", splitblocks=true) # TODO: remove splitbocks=true after warning
+end
+
+function tfim_symbolic_coefficients(coupling_weights::AbstractVector, h::Real)::Vector{Float64}
+  coefficients = Vector{Float64}(undef, length(coupling_weights) + 1)
+  coefficients[1:(end - 1)] .= coupling_weights
+  coefficients[end] = h
+  return coefficients
+end
+
+function test_symbolic_tfim_sample(
+  sym::SymbolicMPO, coupling_weights::AbstractVector, h::Real
+)::Float64
+  
+  symbolic_mpo = instantiate_MPO(
+    sym, tfim_symbolic_coefficients(coupling_weights, h)
+  )
+
+  numeric_mpo = all_to_all_tfim_numeric(sym.sites, coupling_weights, h)
+
+  return norm(add(symbolic_mpo, -numeric_mpo; alg="directsum"))
+end
+
+let N = 50, h = 1.0, num_samples = 5
+  sym = all_to_all_tfim_symbolic(N)
+  num_couplings = num_tfim_couplings(N)
+
+  println(
+    "Constructed a symbolic all-to-all transverse-field Ising MPO template ",
+    "for $N sites",
+  )
+
+  for sample in 1:num_samples
+    error = test_symbolic_tfim_sample(sym, randn(num_couplings), h)
+    println("Sample $sample MPO error: $(round(error; sigdigits=3))")
+    @assert iszero(error)
+  end
+end
+
+nothing
diff --git a/src/MPOConstruction.jl b/src/MPOConstruction.jl
index 6a17302..70ca077 100644
--- a/src/MPOConstruction.jl
+++ b/src/MPOConstruction.jl
@@ -109,6 +109,8 @@ function _remove_symbolic_mpo_boundary_links!(H::MPO, llinks::Vector{<:Index})::
   return nothing
 end
 
+# TODO: Move splitblocks warning to all instantiate_MPO calls
+
 """
     instantiate_MPO(offsets, block_sparse_matrices, sites, llinks; splitblocks, checkflux) -> MPO
 

From 57709d29945fdd2f9db7a1cb21d02ef7d891ce64 Mon Sep 17 00:00:00 2001
From: Ben Corbett <corbett5@llnl.gov>
Date: Mon, 27 Apr 2026 16:25:08 -0600
Subject: [PATCH 11/16] added TFIM example.

---
 examples/symbolic-transverse-field-ising.jl | 9 +++++++--
 1 file changed, 7 insertions(+), 2 deletions(-)

diff --git a/examples/symbolic-transverse-field-ising.jl b/examples/symbolic-transverse-field-ising.jl
index 32110e3..9168b33 100644
--- a/examples/symbolic-transverse-field-ising.jl
+++ b/examples/symbolic-transverse-field-ising.jl
@@ -10,6 +10,7 @@
 
 using ITensors, ITensorMPS, ITensorMPOConstruction
 using Random
+using Test
 
 function num_tfim_couplings(N::Int)::Int
   return (N * (N - 1)) ÷ 2
@@ -89,7 +90,11 @@ function test_symbolic_tfim_sample(
   return norm(add(symbolic_mpo, -numeric_mpo; alg="directsum"))
 end
 
-let N = 50, h = 1.0, num_samples = 5
+@testset "Symbolic TFIM" begin
+  N = 50
+  h = 1.0
+  num_samples = 5
+  
   sym = all_to_all_tfim_symbolic(N)
   num_couplings = num_tfim_couplings(N)
 
@@ -101,7 +106,7 @@ let N = 50, h = 1.0, num_samples = 5
   for sample in 1:num_samples
     error = test_symbolic_tfim_sample(sym, randn(num_couplings), h)
     println("Sample $sample MPO error: $(round(error; sigdigits=3))")
-    @assert iszero(error)
+    @test iszero(error)
   end
 end
 

From 6b4c9419fbb987915ee4bcad331a26a852c7a7fa Mon Sep 17 00:00:00 2001
From: Ben Corbett <corbett5@llnl.gov>
Date: Mon, 27 Apr 2026 23:07:10 -0600
Subject: [PATCH 12/16] Electronic structure example, need to optimize and
 clean.

---
 examples/electronic-structure.jl    | 196 ++++++++++++++++++++++++++--
 src/MPOConstruction.jl              |  37 +++---
 src/SymbolicLocalMatrix.jl          |  42 ++----
 src/large-graph-mpo-common.jl       |  18 +++
 test/symbolic-local-matrix-tests.jl |  13 +-
 5 files changed, 241 insertions(+), 65 deletions(-)

diff --git a/examples/electronic-structure.jl b/examples/electronic-structure.jl
index 876af39..ea11da0 100644
--- a/examples/electronic-structure.jl
+++ b/examples/electronic-structure.jl
@@ -55,12 +55,7 @@ function get_coefficients(N::Int)::Tuple{Array{Float64,2},Array{Float64,4}}
   return h, V
 end
 
-function electronic_structure_OpIDSum(
-  N::Int, h::Array{Float64,2}, V::Array{Float64,4}
-)::Tuple{Vector{<:Index},OpIDSum}
-  ↓ = false
-  ↑ = true
-
+function electronic_structure_sites_and_op_cache(N::Int)
   sites = siteinds("Electron", N; conserve_qns=true)
 
   operatorNames = [
@@ -85,6 +80,17 @@ function electronic_structure_OpIDSum(
   ]
 
   op_cache_vec = to_OpCacheVec(sites, operatorNames)
+  return sites, op_cache_vec
+end
+
+function electronic_structure_OpIDSum(
+  N::Int,
+  h::Array{Float64,2},
+  V::Array{Float64,4},
+  op_cache_vec::OpCacheVec,
+)::OpIDSum
+  ↓ = false
+  ↑ = true
 
   opC(k::Int, spin::Bool) = OpID{UInt8}(2 + spin, k)
   opCdag(k::Int, spin::Bool) = OpID{UInt8}(4 + spin, k)
@@ -129,7 +135,14 @@ function electronic_structure_OpIDSum(
     end
   end
 
-  return sites, os
+  return os
+end
+
+function electronic_structure_OpIDSum(
+  N::Int, h::Array{Float64,2}, V::Array{Float64,4}
+)::Tuple{Vector{<:Index},OpIDSum}
+  sites, op_cache_vec = electronic_structure_sites_and_op_cache(N)
+  return sites, electronic_structure_OpIDSum(N, h, V, op_cache_vec)
 end
 
 for alg in ("VC", "QR")
@@ -162,8 +175,6 @@ for alg in ("VC", "QR")
   end
 end
 
-nothing;
-
 # ## Results
 # Below are the runtime and sparsities for the QR decomposition algorithm and the vertex cover algorithm. Both algorithms produce MPOs of bond dimension ``2 N^2 + 3 N + 2``, which is optimal for a generic set of coefficients. These timings were taken with `julia -t8 --gcthreads=8,1` on a 2021 MacBook Pro with the M1 Max CPU and 32GB of memory.
 #
@@ -192,3 +203,170 @@ nothing;
 # | ``\left[\text{Fe}_2 \text{S} (\text{C H}_3)(\text{SCH}_3)_4 \right]^{3-}`` | 36    | `VC`      | 2702           | 99.26%   |
 #
 # ### Thanks to [Huanchen Zhai](https://scholar.google.com/citations?user=HM_YBL0AAAAJ&hl=en) for providing the discussion and data motivating this section.
+
+# # Symbolic construction
+
+function electronic_structure_fermion_sort_sign(site_ids::Integer...)::Int
+  sign = 1
+  for i in eachindex(site_ids)
+    for j in 1:(i - 1)
+      site_ids[i] < site_ids[j] && (sign = -sign)
+    end
+  end
+  return sign
+end
+
+function electronic_structure_symbolic_OpIDSum(
+  N::Int,
+  op_cache_vec::OpCacheVec,
+)::Tuple{OpIDSum, Vector{Tuple{Int, Int}}, Vector{Tuple{Int, Int, Int, Int, Bool, Bool, Bool, Bool}}}
+  ↓ = false
+  ↑ = true
+
+  opC(k::Int, spin::Bool) = OpID(2 + spin, k)
+  opCdag(k::Int, spin::Bool) = OpID(4 + spin, k)
+
+  os = OpIDSum{4,Int,Int}(2 * N^4 + 2 * N^2, op_cache_vec)
+
+  map_1e = Tuple{Int, Int}[]
+  for p in 1:N
+    for q in 1:N
+      push!(map_1e, (p, q))
+      id = length(map_1e)
+      sign = electronic_structure_fermion_sort_sign(p, q)
+      for spin in (↓, ↑)
+        add!(os, sign * id, opCdag(p, spin), opC(q, spin))
+      end
+    end
+  end
+
+  map_2e = Tuple{Int, Int, Int, Int, Bool, Bool, Bool, Bool}[]
+  for p in 1:N
+    for s_p in (↓, ↑)
+      for q in 1:N, s_q in (↓, ↑)
+        (q, s_q) <= (p, s_p) && continue
+
+        for r in 1:N, s_r in (↓, ↑)
+          for s in 1:N, s_s in (↓, ↑)
+            (s, s_s) <= (r, s_r) && continue
+
+            valid_entry = (s_r == s_q && s_s == s_p) ||
+              (s_s == s_q && s_r == s_p) ||
+              (s_r == s_p && s_s == s_q) ||
+              (s_s == s_p && s_r == s_q)
+
+            !valid_entry && continue
+
+            push!(map_2e, (p, q, r, s, s_p, s_q, s_r, s_s))
+
+            id = length(map_2e) + length(map_1e)
+            sign = electronic_structure_fermion_sort_sign(p, q, r, s)
+            add!(
+              os,
+              sign * id,
+              opCdag(p, s_p),
+              opCdag(q, s_q),
+              opC(r, s_r),
+              opC(s, s_s),
+            )
+          end
+        end
+      end
+    end
+  end
+
+  return os, map_1e, map_2e
+end
+
+function electronic_structure_symbolic_coefficients(
+  h::AbstractMatrix,
+  V::AbstractArray{<:Number,4},
+  map_1e::AbstractVector{<:Tuple{Int,Int}},
+  map_2e::AbstractVector{<:Tuple{Int,Int,Int,Int,Bool,Bool,Bool,Bool}},
+)::Vector{promote_type(eltype(h), eltype(V))}
+  coefficients = Vector{promote_type(eltype(h), eltype(V))}(
+    undef, length(map_1e) + length(map_2e)
+  )
+
+  for (id, (p, q)) in pairs(map_1e)
+    coefficients[id] = electronic_structure_fermion_sort_sign(p, q) * h[p, q]
+  end
+
+  offset = length(map_1e)
+  for (i, (p, q, r, s, s_p, s_q, s_r, s_s)) in pairs(map_2e)
+    coeff = zero(eltype(coefficients))
+    if s_r == s_q && s_s == s_p
+      coeff += V[p, s, q, r]
+    end
+    if s_s == s_q && s_r == s_p
+      coeff -= V[p, r, q, s]
+    end
+    if s_r == s_p && s_s == s_q
+      coeff -= V[q, s, p, r]
+    end
+    if s_s == s_p && s_r == s_q
+      coeff += V[q, r, p, s]
+    end
+
+    coefficients[offset + i] = electronic_structure_fermion_sort_sign(p, q, r, s) * coeff
+  end
+
+  return coefficients
+end
+
+function mpo_relative_difference(A::MPO, B::MPO)::Float64
+  AmB = add(A, -B; alg="directsum")
+  return real(inner(AmB, AmB)) / real(inner(B, B))
+end
+
+let
+  N = 4
+  println("Constructing a symbolic electronic structure MPO for $N sites")
+
+  h, V = get_coefficients(N)
+  sites, op_cache_vec = electronic_structure_sites_and_op_cache(N)
+  numeric_os = electronic_structure_OpIDSum(N, h, V, op_cache_vec)
+  symbolic_os, map_1e, map_2e = electronic_structure_symbolic_OpIDSum(N, op_cache_vec)
+  coefficients = electronic_structure_symbolic_coefficients(h, V, map_1e, map_2e)
+
+  @time "Constructing symbolic MPO" sym = MPO_symbolic(
+    symbolic_os,
+    sites;
+    basis_op_cache_vec=op_cache_vec,
+    check_for_errors=false,
+  )
+  @time "Instantiating symbolic MPO" H_symbolic = instantiate_MPO(
+    sym, coefficients; splitblocks=true, checkflux=false
+  )
+
+  reset_timer!()
+
+  @time "Instantiating symbolic MPO 2" H_symbolic = instantiate_MPO!(
+    H_symbolic, sym, coefficients; checkflux=false
+  )
+
+  print_timer()
+
+  @time "Constructing reference MPO" H_numeric = MPO_new(
+    numeric_os,
+    sites;
+    alg="VC",
+    basis_op_cache_vec=op_cache_vec,
+    splitblocks=true,
+    check_for_errors=false,
+    checkflux=false,
+  )
+
+  relative_difference = mpo_relative_difference(H_symbolic, H_numeric)
+  println("Symbolic instantiation relative difference: $relative_difference")
+  @assert relative_difference < 1e-10
+
+  h2 = 1.1 .* h
+  V2 = 1.1 .* V
+  coefficients2 = electronic_structure_symbolic_coefficients(h2, V2, map_1e, map_2e)
+  H_symbolic2 = instantiate_MPO(sym, coefficients2; splitblocks=true, checkflux=false)
+  @assert mpo_relative_difference(H_symbolic2, H_symbolic) > 0
+  println("Updated symbolic coefficients without rebuilding the symbolic MPO")
+end
+
+nothing;
diff --git a/src/MPOConstruction.jl b/src/MPOConstruction.jl
index 70ca077..85ad849 100644
--- a/src/MPOConstruction.jl
+++ b/src/MPOConstruction.jl
@@ -199,7 +199,7 @@ function _evaluate_symbolic_block_sparse_matrix(
 
   for right_link in eachindex(symbolic_matrix)
     for (left_link, terms) in pairs(symbolic_matrix[right_link])
-      block = convert(Matrix{C}, _evaluate_symbolic_local_matrix(terms, coefficients, op_cache))
+      block = _evaluate_symbolic_local_matrix(C, terms, coefficients, op_cache)
       set!(matrix[right_link], left_link, block)
     end
   end
@@ -241,11 +241,12 @@ function instantiate_MPO(
 
   H = MPO(sym.sites)
 
-  @timeit "to_ITensor" Threads.@threads for n in eachindex(sym.sites)
-    block_sparse_matrices = _evaluate_symbolic_block_sparse_matrices(
+  @timeit "to_ITensor" for n in eachindex(sym.sites)
+    @timeit "_evaluate_symbolic_block_sparse_matrices" block_sparse_matrices = _evaluate_symbolic_block_sparse_matrices(
       sym.block_sparse_matrices[n], coefficients, sym.op_cache_vec[n], C
     )
-    H[n] = to_ITensor(
+
+    @timeit "to_ITensor" H[n] = to_ITensor(
       sym.offsets[n],
       block_sparse_matrices,
       sym.llinks[n],
@@ -287,22 +288,16 @@ function _add_symbolic_mpo_boundary_links!(H::MPO, sym::SymbolicMPO)::Nothing
   return nothing
 end
 
-function _fill_symbolic_mpo_tensor_with_boundary_links!(
-  tensor::ITensor, sym::SymbolicMPO, coefficients::AbstractVector, n::Int
-)::ITensor
-  llink = dag(sym.llinks[n])
-  rlink = sym.llinks[n + 1]
-  site = sym.sites[n]
-  tensor = ITensors.permute(tensor, llink, rlink, prime(site), dag(site); allow_alias=true)
-  tensor .= 0
-
+@timeit function _fill_symbolic_mpo_tensor_with_boundary_links!(
+  tensor, sym::SymbolicMPO, coefficients::AbstractVector, n::Int
+)::Nothing
   op_cache = sym.op_cache_vec[n]
 
   for (offset, matrix) in zip(sym.offsets[n], sym.block_sparse_matrices[n])
     for (right_link, column) in enumerate(matrix)
       shifted_right_link = right_link + offset
       for (left_link, terms) in pairs(column)
-        local_matrix = _evaluate_symbolic_local_matrix(terms, coefficients, op_cache)
+        local_matrix = _evaluate_symbolic_local_matrix(eltype(tensor), terms, coefficients, op_cache)
         for i in axes(local_matrix, 1)
           for j in axes(local_matrix, 2)
             iszero(local_matrix[i, j]) && continue
@@ -313,7 +308,7 @@ function _fill_symbolic_mpo_tensor_with_boundary_links!(
     end
   end
 
-  return tensor
+  return nothing
 end
 
 function _check_symbolic_mpo_template(H::MPO, sym::SymbolicMPO)::Nothing
@@ -361,8 +356,16 @@ function instantiate_MPO!(
 
   _add_symbolic_mpo_boundary_links!(H, sym)
 
-  Threads.@threads for n in eachindex(sym.sites)
-    H[n] = _fill_symbolic_mpo_tensor_with_boundary_links!(H[n], sym, coefficients, n)
+  for n in eachindex(sym.sites)
+    t = H[n]
+    llink = dag(sym.llinks[n])
+    rlink = sym.llinks[n + 1]
+    site = sym.sites[n]
+    t = ITensors.permute(t, llink, rlink, prime(site), dag(site); allow_alias=true)
+    t .= 0
+
+    _fill_symbolic_mpo_tensor_with_boundary_links!(t.tensor, sym, coefficients, n)
+    H[n] = t
   end
 
   _remove_symbolic_mpo_boundary_links!(H, sym.llinks)
diff --git a/src/SymbolicLocalMatrix.jl b/src/SymbolicLocalMatrix.jl
index 580033b..b96f5d1 100644
--- a/src/SymbolicLocalMatrix.jl
+++ b/src/SymbolicLocalMatrix.jl
@@ -149,17 +149,10 @@ function _scale_symbolic_weight(signed_weight_id::Integer, factor)::Int
 end
 
 function _weight_value(signed_weight_id::Integer, coefficients::AbstractVector)
-  _check_symbolic_weight_id(signed_weight_id)
   abs_weight_id = abs(signed_weight_id)
   abs_weight_id == 1 && return signed_weight_id > 0 ? 1.0 : -1.0
 
   user_label = abs_weight_id - 1
-  user_label <= length(coefficients) || throw(
-    ArgumentError(
-      "Missing coefficient for symbolic label $user_label. " *
-      "Received $(length(coefficients)) coefficients.",
-    ),
-  )
 
   value = coefficients[user_label]
   return signed_weight_id > 0 ? value : -value
@@ -176,39 +169,24 @@ function _symbolic_local_matrix_eltype(
 end
 
 function _evaluate_symbolic_local_matrix(
-  terms::SymbolicLocalMatrix, coefficients::AbstractVector, op_cache::Vector{OpInfo}
-)::Matrix
-  isempty(op_cache) && throw(
-    ArgumentError("Cannot evaluate a symbolic local matrix with an empty operator cache."),
-  )
-
-  matrix_size = size(op_cache[1].matrix)
-  result = zeros(_symbolic_local_matrix_eltype(coefficients, op_cache), matrix_size)
+  ::Type{C}, terms::SymbolicLocalMatrix, coefficients::AbstractVector, op_cache::Vector{OpInfo}
+)::Matrix{C} where {C}
+  result = zeros(C, size(op_cache[1].matrix))
 
+  needs_jw = 0
   for (signed_weight_id, signed_local_op_id) in terms
-    _check_symbolic_weight_id(signed_weight_id)
-    _check_symbolic_local_op_id(signed_local_op_id)
-
     local_op_id = abs(Int(signed_local_op_id))
-    local_op_id <= length(op_cache) || throw(
-      ArgumentError(
-        "Symbolic local operator id $local_op_id is not available in a cache with " *
-        "$(length(op_cache)) operators.",
-      ),
-    )
-
     local_op = op_cache[local_op_id].matrix
-    size(local_op) == matrix_size || throw(
-      ArgumentError(
-        "Operators in a symbolic local matrix must have matching matrix sizes. " *
-        "Expected $matrix_size, got $(size(local_op)).",
-      ),
-    )
+    weight = _weight_value(signed_weight_id, coefficients)
 
+    needs_jw += signed_local_op_id < 0
     add_to_local_matrix!(
-      result, _weight_value(signed_weight_id, coefficients), local_op, signed_local_op_id < 0
+      result, weight, local_op, false
     )
   end
 
+  @assert needs_jw ∈ (0, length(terms))
+  needs_jw > 0 && apply_jw_string!(result)
+
   return result
 end
diff --git a/src/large-graph-mpo-common.jl b/src/large-graph-mpo-common.jl
index 10b9a4c..bb152a1 100644
--- a/src/large-graph-mpo-common.jl
+++ b/src/large-graph-mpo-common.jl
@@ -285,6 +285,24 @@ carry the remaining operator tuples.
   return g
 end
 
+# TODO: Replace add_to_local_matrix! with this
+function apply_jw_string!(a::Matrix)::Nothing
+  if size(a, 1) == 2
+    @inbounds for i in 1:2
+      a[i, 2] *= -1
+    end
+  elseif size(a, 1) == 4
+    @inbounds for i in 1:4
+      a[i, 2] *= -1
+      a[i, 3] *= -1
+    end
+  else
+    error("Unknown fermionic site.")
+  end
+
+  return nothing
+end
+
 """
     add_to_local_matrix!(a, weight, local_op, needs_JW_string) -> Nothing
 
diff --git a/test/symbolic-local-matrix-tests.jl b/test/symbolic-local-matrix-tests.jl
index 21b66bc..72064ca 100644
--- a/test/symbolic-local-matrix-tests.jl
+++ b/test/symbolic-local-matrix-tests.jl
@@ -26,7 +26,6 @@ function test_symbolic_local_matrix_operations()::Nothing
   @test _weight_value(-1, coefficients) == -1.0
   @test _weight_value(3, coefficients) == coefficients[2]
   @test _weight_value(-3, coefficients) == -coefficients[2]
-  @test_throws ArgumentError _weight_value(4, coefficients)
 
   @test _internal_symbolic_id(1) == 2
   @test _internal_symbolic_id(3) == 4
@@ -60,16 +59,16 @@ function test_symbolic_local_matrix_operations()::Nothing
   qubit_op_cache_vec = to_OpCacheVec(qubit_sites, [["I", "X"]])
   duplicate_terms = SymbolicLocalMatrix{Int}([(3, 2), (3, 2)])
   @test _evaluate_symbolic_local_matrix(
-    duplicate_terms, coefficients, qubit_op_cache_vec[1]
+    Float64, duplicate_terms, coefficients, qubit_op_cache_vec[1]
   ) ≈ 2 * coefficients[2] * qubit_op_cache_vec[1][2].matrix
 
   fermion_sites = siteinds("Fermion", 1)
   fermion_op_cache_vec = to_OpCacheVec(fermion_sites, [["I", "C", "Cdag", "N"]])
   plain_matrix = _evaluate_symbolic_local_matrix(
-    SymbolicLocalMatrix{Int}([(3, 2)]), coefficients, fermion_op_cache_vec[1]
+    Float64, SymbolicLocalMatrix{Int}([(3, 2)]), coefficients, fermion_op_cache_vec[1]
   )
   jw_matrix = _evaluate_symbolic_local_matrix(
-    SymbolicLocalMatrix{Int}([(3, -2)]), coefficients, fermion_op_cache_vec[1]
+    Float64, SymbolicLocalMatrix{Int}([(3, -2)]), coefficients, fermion_op_cache_vec[1]
   )
 
   expected_plain_matrix = coefficients[2] * fermion_op_cache_vec[1][2].matrix
@@ -274,7 +273,7 @@ function test_symbolic_mpo_graph_preserves_duplicate_signed_ids()::Nothing
   @test terms == [(2, 2), (3, 2)]
 
   coefficients = [11.0, 17.0]
-  evaluated = _evaluate_symbolic_local_matrix(terms, coefficients, op_cache_vec[1])
+  evaluated = _evaluate_symbolic_local_matrix(Float64, terms, coefficients, op_cache_vec[1])
   @test evaluated ≈ sum(coefficients) * op_cache_vec[1][2].matrix
 
   return nothing
@@ -296,7 +295,7 @@ function test_symbolic_mpo_graph_preserves_duplicate_label_multiplicity()::Nothi
   @test terms == [(2, 2), (2, 2)]
 
   coefficient = 13.0
-  evaluated = _evaluate_symbolic_local_matrix(terms, [coefficient], op_cache_vec[1])
+  evaluated = _evaluate_symbolic_local_matrix(Float64, terms, [coefficient], op_cache_vec[1])
   @test evaluated ≈ 2 * coefficient * op_cache_vec[1][2].matrix
 
   return nothing
@@ -339,7 +338,7 @@ function test_symbolic_mpo_graph_cancels_opposite_signed_ids()::Nothing
   terms = symbolic_terms_from_first_left_vertex(g)
   @test isempty(terms)
 
-  evaluated = _evaluate_symbolic_local_matrix(terms, [7.0], basis_op_cache_vec[1])
+  evaluated = _evaluate_symbolic_local_matrix(Float64, terms, [7.0], basis_op_cache_vec[1])
   @test iszero(evaluated)
 
   return nothing

From 6287a2e035207acbaab4c34b1e9176dc1f87b340 Mon Sep 17 00:00:00 2001
From: Ben Corbett <corbett5@llnl.gov>
Date: Tue, 28 Apr 2026 10:02:43 -0600
Subject: [PATCH 13/16] Removed TODO

---
 src/large-graph-mpo-vertex-cover.jl | 2 --
 1 file changed, 2 deletions(-)

diff --git a/src/large-graph-mpo-vertex-cover.jl b/src/large-graph-mpo-vertex-cover.jl
index 65b55e2..67d2dc3 100644
--- a/src/large-graph-mpo-vertex-cover.jl
+++ b/src/large-graph-mpo-vertex-cover.jl
@@ -121,7 +121,6 @@ otherwise this also builds `next_edges_of_cc` for the graph at site `n + 1`.
     end
   end
 
-  # TODO: Consider nested multithreading
   Threads.@threads for cc in 1:num_connected_components(ccs)
     lvs_of_component::Vector{Int} = ccs.lvs_of_component[cc]
     position_of_rvs_in_component = ccs.position_of_rvs_in_component
@@ -145,7 +144,6 @@ otherwise this also builds `next_edges_of_cc` for the graph at site `n + 1`.
     end
 
     ## Construct the tensor from the right cover.
-    # TODO: Merge this loop with the one above, using `in_left_cover`
     let
       in_left_cover = falses(length(lvs_of_component))
       @inbounds for local_id in left_cover

From 4865b5f094577259a2caf74ee5fefa460b503289 Mon Sep 17 00:00:00 2001
From: Ben Corbett <corbett5@llnl.gov>
Date: Tue, 28 Apr 2026 11:02:04 -0600
Subject: [PATCH 14/16] Faster, but not as fast as it should be and messy.

---
 examples/electronic-structure.jl | 109 ++++++++-------
 src/MPOConstruction.jl           | 231 +++++++++++++++++++++++++++++--
 src/SymbolicLocalMatrix.jl       |  19 ++-
 test/symbolic-mpo-tests.jl       |   9 ++
 4 files changed, 299 insertions(+), 69 deletions(-)

diff --git a/examples/electronic-structure.jl b/examples/electronic-structure.jl
index ea11da0..72cba6b 100644
--- a/examples/electronic-structure.jl
+++ b/examples/electronic-structure.jl
@@ -145,35 +145,35 @@ function electronic_structure_OpIDSum(
   return sites, electronic_structure_OpIDSum(N, h, V, op_cache_vec)
 end
 
-for alg in ("VC", "QR")
-  for N in [10, 20]
-    println("Constructing the electronic structure MPO for $N sites using $alg")
-
-    reset_timer!()
-    @time "Constructing OpIDSum" sites, os = electronic_structure_OpIDSum(
-      N, get_coefficients(N)...
-    )
-    @time "Constructing MPO" H = MPO_new(
-      os,
-      sites;
-      alg,
-      basis_op_cache_vec=os.op_cache_vec,
-      splitblocks=true,
-      check_for_errors=false,
-      checkflux=false,
-    ) # TODO: remove splitblocks=true after warning
-    N > 5 && print_timer()
-
-    percent_sparse = round(100 * sparsity(H); digits=2)
-    println(
-      "The maximum bond dimension is $(maxlinkdim(H)), sparsity = $percent_sparse%",
-    )
-    @assert maxlinkdim(H) == 2 * N^2 + 3 * N + 2
-
-    GC.gc(true)
-    println()
-  end
-end
+# for alg in ("VC", "QR")
+#   for N in [10, 20]
+#     println("Constructing the electronic structure MPO for $N sites using $alg")
+
+#     reset_timer!()
+#     @time "Constructing OpIDSum" sites, os = electronic_structure_OpIDSum(
+#       N, get_coefficients(N)...
+#     )
+#     @time "Constructing MPO" H = MPO_new(
+#       os,
+#       sites;
+#       alg,
+#       basis_op_cache_vec=os.op_cache_vec,
+#       splitblocks=true,
+#       check_for_errors=false,
+#       checkflux=false,
+#     ) # TODO: remove splitblocks=true after warning
+#     N > 5 && print_timer()
+
+#     percent_sparse = round(100 * sparsity(H); digits=2)
+#     println(
+#       "The maximum bond dimension is $(maxlinkdim(H)), sparsity = $percent_sparse%",
+#     )
+#     @assert maxlinkdim(H) == 2 * N^2 + 3 * N + 2
+
+#     GC.gc(true)
+#     println()
+#   end
+# end
 
 # ## Results
 # Below are the runtime and sparsities for the QR decomposition algorithm and the vertex cover algorithm. Both algorithms produce MPOs of bond dimension ``2 N^2 + 3 N + 2``, which is optimal for a generic set of coefficients. These timings were taken with `julia -t8 --gcthreads=8,1` on a 2021 MacBook Pro with the M1 Max CPU and 32GB of memory.
@@ -320,7 +320,7 @@ function mpo_relative_difference(A::MPO, B::MPO)::Float64
 end
 
 let
-  N = 4
+  N = 50
   println("Constructing a symbolic electronic structure MPO for $N sites")
 
   h, V = get_coefficients(N)
@@ -335,10 +335,19 @@ let
     basis_op_cache_vec=op_cache_vec,
     check_for_errors=false,
   )
+
+  H_symbolic = instantiate_MPO(
+    sym, coefficients; splitblocks=true, checkflux=false
+  )
+
   @time "Instantiating symbolic MPO" H_symbolic = instantiate_MPO(
     sym, coefficients; splitblocks=true, checkflux=false
   )
 
+  instantiate_MPO!(
+    H_symbolic, sym, coefficients; checkflux=false
+  )
+
   reset_timer!()
 
   @time "Instantiating symbolic MPO 2" H_symbolic = instantiate_MPO!(
@@ -347,26 +356,26 @@ let
 
   print_timer()
 
-  @time "Constructing reference MPO" H_numeric = MPO_new(
-    numeric_os,
-    sites;
-    alg="VC",
-    basis_op_cache_vec=op_cache_vec,
-    splitblocks=true,
-    check_for_errors=false,
-    checkflux=false,
-  )
-
-  relative_difference = mpo_relative_difference(H_symbolic, H_numeric)
-  println("Symbolic instantiation relative difference: $relative_difference")
-  @assert relative_difference < 1e-10
-
-  h2 = 1.1 .* h
-  V2 = 1.1 .* V
-  coefficients2 = electronic_structure_symbolic_coefficients(h2, V2, map_1e, map_2e)
-  H_symbolic2 = instantiate_MPO(sym, coefficients2; splitblocks=true, checkflux=false)
-  @assert mpo_relative_difference(H_symbolic2, H_symbolic) > 0
-  println("Updated symbolic coefficients without rebuilding the symbolic MPO")
+  # @time "Constructing reference MPO" H_numeric = MPO_new(
+  #   numeric_os,
+  #   sites;
+  #   alg="VC",
+  #   basis_op_cache_vec=op_cache_vec,
+  #   splitblocks=true,
+  #   check_for_errors=false,
+  #   checkflux=false,
+  # )
+
+  # relative_difference = mpo_relative_difference(H_symbolic, H_numeric)
+  # println("Symbolic instantiation relative difference: $relative_difference")
+  # @assert relative_difference < 1e-10
+
+  # h2 = 1.1 .* h
+  # V2 = 1.1 .* V
+  # coefficients2 = electronic_structure_symbolic_coefficients(h2, V2, map_1e, map_2e)
+  # H_symbolic2 = instantiate_MPO(sym, coefficients2; splitblocks=true, checkflux=false)
+  # @assert mpo_relative_difference(H_symbolic2, H_symbolic) > 0
+  # println("Updated symbolic coefficients without rebuilding the symbolic MPO")
 end
 
 nothing;
diff --git a/src/MPOConstruction.jl b/src/MPOConstruction.jl
index 85ad849..83e41e7 100644
--- a/src/MPOConstruction.jl
+++ b/src/MPOConstruction.jl
@@ -131,7 +131,7 @@ function instantiate_MPO(
 )::MPO where {ValType<:Number}
   H = MPO(sites)
 
-  @timeit "to_ITensor" Threads.@threads for n in 1:length(sites)
+  @timeit "to_ITensor" for n in 1:length(sites)
     H[n] = to_ITensor(
       offsets[n], block_sparse_matrices[n], llinks[n], llinks[n + 1], sites[n]; splitblocks
     )
@@ -288,20 +288,220 @@ function _add_symbolic_mpo_boundary_links!(H::MPO, sym::SymbolicMPO)::Nothing
   return nothing
 end
 
-@timeit function _fill_symbolic_mpo_tensor_with_boundary_links!(
-  tensor, sym::SymbolicMPO, coefficients::AbstractVector, n::Int
-)::Nothing
+function _fill_symbolic_mpo_template_tensor!(
+  tensor::NDTensors.DenseTensor{C,4},
+  sym::SymbolicMPO,
+  coefficients::AbstractVector,
+  n::Int,
+)::Nothing where {C}
+  op_cache = sym.op_cache_vec[n]
+  tensor_data = array(tensor)
+
+  fill!(tensor_data, zero(C))
+  @timeit "evaluate and fill entries" for (offset, matrix) in zip(
+    sym.offsets[n], sym.block_sparse_matrices[n]
+  )
+    for (right_link, column) in enumerate(matrix)
+      shifted_right_link = right_link + offset
+      for (left_link, terms) in pairs(column)
+        local_matrix = _evaluate_symbolic_local_matrix(C, terms, coefficients, op_cache)
+        @inbounds for site_prime_link in axes(local_matrix, 1),
+          site_link in axes(local_matrix, 2)
+
+          value = local_matrix[site_prime_link, site_link]
+          iszero(value) && continue
+          tensor_data[left_link, shifted_right_link, site_prime_link, site_link] = value
+        end
+      end
+    end
+  end
+
+  return nothing
+end
+
+@inline function _block_key(block::Block{4})::NTuple{4,Int}
+  return (Int(block[1]), Int(block[2]), Int(block[3]), Int(block[4]))
+end
+
+@inline function _block_storage_index(
+  block_offset::Int, entry_offsets::NTuple{4,Int}, block_dims::NTuple{4,Int}
+)::Int
+  storage_index = block_offset + entry_offsets[1]
+  stride = block_dims[1]
+  for dim in 2:4
+    storage_index += (entry_offsets[dim] - 1) * stride
+    stride *= block_dims[dim]
+  end
+  return storage_index
+end
+
+function _try_fill_symbolic_mpo_splitblocks_sequential!(
+  tensor::NDTensors.BlockSparseTensor{C,4},
+  sym::SymbolicMPO,
+  coefficients::AbstractVector,
+  n::Int,
+  left_blocks::Vector{Int},
+  right_blocks::Vector{Int},
+  site_prime_blocks::Vector{Int},
+  site_blocks::Vector{Int},
+  tensor_data,
+)::Bool where {C}
+  length(tensor_data) == NDTensors.nnzblocks(tensor) || return false
+
+  op_cache = sym.op_cache_vec[n]
+  block_iterator = ITensors.eachnzblock(tensor)
+  block_state = iterate(block_iterator)
+  current_block_key = isnothing(block_state) ? nothing : _block_key(block_state[1])
+  storage_index = 1
+
+  fill!(tensor_data, zero(C))
+
+  local_matrix = zeros(C, size(op_cache[1].matrix))
+  @timeit "sequential splitblock refill" for (offset, matrix) in zip(
+    sym.offsets[n], sym.block_sparse_matrices[n]
+  )
+    for (right_link, column) in enumerate(matrix)
+      shifted_right_link = right_link + offset
+      right_block = right_blocks[shifted_right_link]
+      for (left_link, terms) in pairs(column)
+        left_block = left_blocks[left_link]
+        _evaluate_symbolic_local_matrix!(local_matrix, terms, coefficients, op_cache)
+
+        @inbounds for site_prime_link in axes(local_matrix, 1),
+          site_link in axes(local_matrix, 2)
+
+          expected_block_key = (
+            left_block,
+            right_block,
+            site_prime_blocks[site_prime_link],
+            site_blocks[site_link],
+          )
+          value = local_matrix[site_prime_link, site_link]
+
+          if !isnothing(current_block_key) && expected_block_key == current_block_key
+            tensor_data[storage_index] = value
+            storage_index += 1
+            block_state = iterate(block_iterator, block_state[2])
+            current_block_key = isnothing(block_state) ? nothing : _block_key(block_state[1])
+          elseif !iszero(value)
+            return false
+          end
+        end
+      end
+    end
+  end
+
+  return isnothing(current_block_key)
+end
+
+function _fill_symbolic_mpo_template_tensor!(
+  tensor::NDTensors.BlockSparseTensor{C,4},
+  sym::SymbolicMPO,
+  coefficients::AbstractVector,
+  n::Int,
+)::Nothing where {C}
+  tensor_inds = inds(tensor)
+
+  left_blocks, left_offsets = _qn_position_map(tensor_inds[1])
+  right_blocks, right_offsets = _qn_position_map(tensor_inds[2])
+  site_prime_blocks, site_prime_offsets = _qn_position_map(tensor_inds[3])
+  site_blocks, site_offsets = _qn_position_map(tensor_inds[4])
+
+  tensor_data = NDTensors.data(tensor)
+  if _try_fill_symbolic_mpo_splitblocks_sequential!(
+    tensor,
+    sym,
+    coefficients,
+    n,
+    left_blocks,
+    right_blocks,
+    site_prime_blocks,
+    site_blocks,
+    tensor_data,
+  )
+    return nothing
+  end
+
+  block_offsets = sizehint!(Dict{NTuple{4,Int},Int}(), NDTensors.nnzblocks(tensor))
+
+  
+  fill!(tensor_data, zero(C))
+
+  @timeit "index template blocks" begin
+    for block in ITensors.eachnzblock(tensor)
+      block_offsets[_block_key(block)] = NDTensors.offset(NDTensors.blockoffsets(tensor), block)
+    end
+  end
+
+  @timeit "evaluate and fill entries" _fill_symbolic_mpo_tensor_entries!(
+    tensor,
+    sym,
+    coefficients,
+    n,
+    left_blocks,
+    right_blocks,
+    site_prime_blocks,
+    site_blocks,
+    left_offsets,
+    right_offsets,
+    site_prime_offsets,
+    site_offsets,
+    tensor_inds,
+    tensor_data,
+    block_offsets,
+  )
+
+  return nothing
+end
+
+function _fill_symbolic_mpo_tensor_entries!(
+  tensor::NDTensors.BlockSparseTensor{C,4},
+  sym::SymbolicMPO,
+  coefficients::AbstractVector,
+  n::Int,
+  left_blocks::Vector{Int},
+  right_blocks::Vector{Int},
+  site_prime_blocks::Vector{Int},
+  site_blocks::Vector{Int},
+  left_offsets::Vector{Int},
+  right_offsets::Vector{Int},
+  site_prime_offsets::Vector{Int},
+  site_offsets::Vector{Int},
+  tensor_inds,
+  tensor_data,
+  block_offsets::Dict{NTuple{4,Int},Int},
+)::Nothing where {C}
   op_cache = sym.op_cache_vec[n]
 
   for (offset, matrix) in zip(sym.offsets[n], sym.block_sparse_matrices[n])
     for (right_link, column) in enumerate(matrix)
       shifted_right_link = right_link + offset
+      right_block = right_blocks[shifted_right_link]
+      right_offset = right_offsets[shifted_right_link]
       for (left_link, terms) in pairs(column)
-        local_matrix = _evaluate_symbolic_local_matrix(eltype(tensor), terms, coefficients, op_cache)
-        for i in axes(local_matrix, 1)
-          for j in axes(local_matrix, 2)
-            iszero(local_matrix[i, j]) && continue
-            tensor[left_link, shifted_right_link, i, j] = local_matrix[i, j]
+        left_block = left_blocks[left_link]
+        left_offset = left_offsets[left_link]
+        local_matrix = _evaluate_symbolic_local_matrix(C, terms, coefficients, op_cache)
+
+        @inbounds for site_prime_link in axes(local_matrix, 1),
+          site_link in axes(local_matrix, 2)
+
+          value = local_matrix[site_prime_link, site_link]
+          iszero(value) && continue
+
+          site_prime_block = site_prime_blocks[site_prime_link]
+          site_block = site_blocks[site_link]
+          block_id = (left_block, right_block, site_prime_block, site_block)
+
+          block_offset = get(block_offsets, block_id, nothing)
+          if isnothing(block_offset)
+            tensor[left_link, shifted_right_link, site_prime_link, site_link] = value
+          else
+            site_prime_offset = site_prime_offsets[site_prime_link]
+            site_offset = site_offsets[site_link]
+            entry_offsets = (left_offset, right_offset, site_prime_offset, site_offset)
+            block_dims = ntuple(i -> blockdim(tensor_inds[i], block_id[i]), Val(4))
+            tensor_data[_block_storage_index(block_offset, entry_offsets, block_dims)] = value
           end
         end
       end
@@ -356,19 +556,22 @@ function instantiate_MPO!(
 
   _add_symbolic_mpo_boundary_links!(H, sym)
 
-  for n in eachindex(sym.sites)
+  @timeit "refill template tensors" for n in eachindex(sym.sites)
     t = H[n]
     llink = dag(sym.llinks[n])
     rlink = sym.llinks[n + 1]
     site = sym.sites[n]
-    t = ITensors.permute(t, llink, rlink, prime(site), dag(site); allow_alias=true)
-    t .= 0
+    t = ITensors.permute(
+      t, llink, rlink, prime(site), dag(site); allow_alias=true
+    )
 
-    _fill_symbolic_mpo_tensor_with_boundary_links!(t.tensor, sym, coefficients, n)
+    _fill_symbolic_mpo_template_tensor!(t.tensor, sym, coefficients, n)
     H[n] = t
   end
 
-  _remove_symbolic_mpo_boundary_links!(H, sym.llinks)
+  _remove_symbolic_mpo_boundary_links!(
+    H, sym.llinks
+  )
 
   if checkflux
     @timeit "checkflux" Threads.@threads for n in eachindex(sym.sites)
diff --git a/src/SymbolicLocalMatrix.jl b/src/SymbolicLocalMatrix.jl
index b96f5d1..cf29f7c 100644
--- a/src/SymbolicLocalMatrix.jl
+++ b/src/SymbolicLocalMatrix.jl
@@ -168,10 +168,10 @@ function _symbolic_local_matrix_eltype(
   return val_type
 end
 
-function _evaluate_symbolic_local_matrix(
-  ::Type{C}, terms::SymbolicLocalMatrix, coefficients::AbstractVector, op_cache::Vector{OpInfo}
-)::Matrix{C} where {C}
-  result = zeros(C, size(op_cache[1].matrix))
+function _evaluate_symbolic_local_matrix!(
+  result::Matrix{C}, terms::SymbolicLocalMatrix, coefficients::AbstractVector, op_cache::Vector{OpInfo}
+)::Nothing where {C}
+  result .= 0
 
   needs_jw = 0
   for (signed_weight_id, signed_local_op_id) in terms
@@ -188,5 +188,14 @@ function _evaluate_symbolic_local_matrix(
   @assert needs_jw ∈ (0, length(terms))
   needs_jw > 0 && apply_jw_string!(result)
 
-  return result
+  return nothing
 end
+
+# TODO: Make sure this isn't used anywhere
+function _evaluate_symbolic_local_matrix(
+  ::Type{C}, terms::SymbolicLocalMatrix, coefficients::AbstractVector, op_cache::Vector{OpInfo}
+)::Matrix{C} where {C}
+  result = Matrix{C}(undef, size(op_cache[1].matrix))
+  _evaluate_symbolic_local_matrix!(result, terms, coefficients, op_cache)
+  return result
+end
\ No newline at end of file
diff --git a/test/symbolic-mpo-tests.jl b/test/symbolic-mpo-tests.jl
index a4c5625..6a106e8 100644
--- a/test/symbolic-mpo-tests.jl
+++ b/test/symbolic-mpo-tests.jl
@@ -228,6 +228,15 @@ function test_symbolic_mpo_inplace_instantiation()::Nothing
   @test linkdims(H) == original_linkdims
   test_mpos_approx_equal(H, expected)
 
+  coefficients3 = [0.0, 1.4, 0.6]
+  result = instantiate_MPO!(H, sym, coefficients3; checkflux=true)
+  expected = instantiate_MPO(sym, coefficients3; splitblocks=true, checkflux=true)
+
+  @test result === H
+  @test objectid(H) == original_objectid
+  @test linkdims(H) == original_linkdims
+  test_mpos_approx_equal(H, expected)
+
   @test_throws MethodError instantiate_MPO!(H, sym, coefficients2; splitblocks=false)
 
   return nothing

From f837d611a1ca7099618434bfa1ba8732c9a6b31f Mon Sep 17 00:00:00 2001
From: Ben Corbett <corbett5@llnl.gov>
Date: Wed, 29 Apr 2026 14:11:43 -0600
Subject: [PATCH 15/16] Unified BlockSparseMatrix

---
 examples/fermi-hubbard-tc.jl        |  4 +++-
 src/MPOConstruction.jl              | 22 +++++++++++++---------
 src/SymbolicLocalMatrix.jl          | 11 +----------
 src/large-graph-mpo-common.jl       | 16 ++++++++++------
 src/large-graph-mpo-qr.jl           |  4 ++--
 src/large-graph-mpo-vertex-cover.jl | 20 ++++++++++----------
 src/sparse_tensor_construction.jl   | 24 ++++++++++++------------
 test/tensor-construction-tests.jl   |  2 +-
 8 files changed, 52 insertions(+), 51 deletions(-)

diff --git a/examples/fermi-hubbard-tc.jl b/examples/fermi-hubbard-tc.jl
index 91d9f0e..ad51e83 100644
--- a/examples/fermi-hubbard-tc.jl
+++ b/examples/fermi-hubbard-tc.jl
@@ -431,7 +431,9 @@ using Serialization
 function call_back(
   n::Int,
   offsets::Vector{Vector{Int}},
-  block_sparse_matrices::Vector{Vector{ITensorMPOConstruction.BlockSparseMatrix{ValType}}},
+  block_sparse_matrices::Vector{
+    Vector{ITensorMPOConstruction.BlockSparseMatrix{Matrix{ValType}}}
+  },
   sites::Vector{<:Index},
   llinks::Vector{<:Index},
   g::ITensorMPOConstruction.MPOGraph,
diff --git a/src/MPOConstruction.jl b/src/MPOConstruction.jl
index 83e41e7..cff29e8 100644
--- a/src/MPOConstruction.jl
+++ b/src/MPOConstruction.jl
@@ -123,7 +123,7 @@ Afterward, the dummy left and right boundary links are contracted away. If
 """
 function instantiate_MPO(
   offsets::Vector{Vector{Int}},
-  block_sparse_matrices::Vector{Vector{BlockSparseMatrix{ValType}}},
+  block_sparse_matrices::Vector{Vector{BlockSparseMatrix{Matrix{ValType}}}},
   sites::Vector{<:Index},
   llinks::Vector{<:Index};
   splitblocks::Bool,
@@ -162,7 +162,7 @@ struct SymbolicMPO{
   Ti<:Integer,Sites<:AbstractVector{<:Index},Links<:AbstractVector{<:Index}
 }
   offsets::Vector{Vector{Int}}
-  block_sparse_matrices::Vector{Vector{SymbolicBlockSparseMatrix{Ti}}}
+  block_sparse_matrices::Vector{Vector{BlockSparseMatrix{SymbolicLocalMatrix{Ti}}}}
   sites::Sites
   llinks::Links
   op_cache_vec::OpCacheVec
@@ -190,11 +190,11 @@ function _symbolic_mpo_eltype(sym::SymbolicMPO, coefficients::AbstractVector)::T
 end
 
 function _evaluate_symbolic_block_sparse_matrix(
-  symbolic_matrix::SymbolicBlockSparseMatrix{Ti},
+  symbolic_matrix::BlockSparseMatrix{SymbolicLocalMatrix{Ti}},
   coefficients::AbstractVector,
   op_cache::Vector{OpInfo},
   ::Type{C},
-)::BlockSparseMatrix{C} where {Ti,C}
+)::BlockSparseMatrix{Matrix{C}} where {Ti,C}
   matrix = [Dictionary{Int,Matrix{C}}() for _ in eachindex(symbolic_matrix)]
 
   for right_link in eachindex(symbolic_matrix)
@@ -208,12 +208,12 @@ function _evaluate_symbolic_block_sparse_matrix(
 end
 
 function _evaluate_symbolic_block_sparse_matrices(
-  symbolic_matrices::Vector{SymbolicBlockSparseMatrix{Ti}},
+  symbolic_matrices::Vector{BlockSparseMatrix{SymbolicLocalMatrix{Ti}}},
   coefficients::AbstractVector,
   op_cache::Vector{OpInfo},
   ::Type{C},
-)::Vector{BlockSparseMatrix{C}} where {Ti,C}
-  matrices = Vector{BlockSparseMatrix{C}}(undef, length(symbolic_matrices))
+)::Vector{BlockSparseMatrix{Matrix{C}}} where {Ti,C}
+  matrices = Vector{BlockSparseMatrix{Matrix{C}}}(undef, length(symbolic_matrices))
   for cc in eachindex(symbolic_matrices)
     matrices[cc] = _evaluate_symbolic_block_sparse_matrix(
       symbolic_matrices[cc], coefficients, op_cache, C
@@ -654,7 +654,9 @@ function MPO_symbolic(
   end
 
   offsets = Vector{Vector{Int}}(undef, length(sites))
-  block_sparse_matrices = Vector{Vector{SymbolicBlockSparseMatrix{Ti}}}(undef, length(sites))
+  block_sparse_matrices = Vector{Vector{BlockSparseMatrix{SymbolicLocalMatrix{Ti}}}}(
+    undef, length(sites)
+  )
 
   @time_if output_level 0 "Constructing symbolic MPO terms" resume_MPO_construction!(
     1,
@@ -734,7 +736,9 @@ function MPO_new(
   end
 
   offsets = Vector{Vector{Int}}(undef, length(sites))
-  block_sparse_matrices = Vector{Vector{BlockSparseMatrix{ValType}}}(undef, length(sites))
+  block_sparse_matrices = Vector{Vector{BlockSparseMatrix{Matrix{ValType}}}}(
+    undef, length(sites)
+  )
 
   @time_if output_level 0 "Constructing MPO terms" resume_MPO_construction!(
     1,
diff --git a/src/SymbolicLocalMatrix.jl b/src/SymbolicLocalMatrix.jl
index cf29f7c..150cca4 100644
--- a/src/SymbolicLocalMatrix.jl
+++ b/src/SymbolicLocalMatrix.jl
@@ -12,15 +12,6 @@ should be applied while evaluating the cached operator matrix.
 """
 const SymbolicLocalMatrix{Ti<:Integer} = Vector{Tuple{Int,Ti}}
 
-"""
-    SymbolicBlockSparseMatrix{Ti}
-
-Block-sparse symbolic storage mirroring `BlockSparseMatrix`, but with each
-block entry represented as a [`SymbolicLocalMatrix`](@ref).
-"""
-const SymbolicBlockSparseMatrix{Ti<:Integer} =
-  Vector{Dictionary{Int,SymbolicLocalMatrix{Ti}}}
-
 function _check_symbolic_weight_id(signed_weight_id::Integer)::Nothing
   signed_weight_id == 0 && throw(
     ArgumentError(
@@ -198,4 +189,4 @@ function _evaluate_symbolic_local_matrix(
   result = Matrix{C}(undef, size(op_cache[1].matrix))
   _evaluate_symbolic_local_matrix!(result, terms, coefficients, op_cache)
   return result
-end
\ No newline at end of file
+end
diff --git a/src/large-graph-mpo-common.jl b/src/large-graph-mpo-common.jl
index bb152a1..664baf9 100644
--- a/src/large-graph-mpo-common.jl
+++ b/src/large-graph-mpo-common.jl
@@ -1,15 +1,17 @@
 """
-    BlockSparseMatrix{C}
+    BlockSparseMatrix{MatrixType}
 
 Vector-backed block-sparse matrix representation used for intermediate MPO
 tensor storage.
 
 The outer vector is indexed by component-local `right_link`. Each inner
-`Dictionary` maps `left_link` to the dense local operator matrix for that block.
+`Dictionary` maps `left_link` to one local block value. Numeric construction
+uses dense local operator matrices as `BlockSparseMatrix{Matrix{C}}`, while the
+symbolic path stores term lists as `BlockSparseMatrix{SymbolicLocalMatrix{Ti}}`.
 `at_site!` later returns offsets that place component-local right-link ids into
 the full outgoing MPO bond.
 """
-BlockSparseMatrix{C} = Vector{Dictionary{Int,Matrix{C}}}
+BlockSparseMatrix{MatrixType} = Vector{Dictionary{Int,MatrixType}}
 
 """
     MPOGraph{N,C,Ti}
@@ -367,7 +369,9 @@ function merge_qn_sectors(
   return new_order, new_qi
 end
 
-function _check_qr_block_storage(::Type{SymbolicBlockSparseMatrix{Ti}})::Nothing where {Ti}
+function _check_qr_block_storage(
+  ::Type{BlockSparseMatrix{SymbolicLocalMatrix{Ti}}}
+)::Nothing where {Ti}
   throw(ArgumentError("QR construction is only supported for numeric block storage."))
 end
 
@@ -376,9 +380,9 @@ function _check_qr_block_storage(::Type)::Nothing
 end
 
 """
-    at_site!(ValType, g, n, sites, tol, absolute_tol, op_cache_vec, alg;
+    at_site!(MatrixType, g, n, sites, tol, absolute_tol, op_cache_vec, alg;
              combine_qn_sectors, output_level=0)
-        -> Tuple{MPOGraph,Vector{Int},Vector{BlockSparseMatrix{ValType}},Index}
+        -> Tuple{MPOGraph,Vector{Int},Vector{MatrixType},Index}
 
 Process one site of the MPO construction algorithm.
 
diff --git a/src/large-graph-mpo-qr.jl b/src/large-graph-mpo-qr.jl
index db2e1ac..edef29a 100644
--- a/src/large-graph-mpo-qr.jl
+++ b/src/large-graph-mpo-qr.jl
@@ -93,7 +93,7 @@ outgoing edges for the next site. For nonterminal sites, the consumed left
 adjacency list is cleared after the outgoing edges are built.
 """
 function process_single_left_vertex_cc!(
-  matrix_of_cc::Vector{BlockSparseMatrix{ValType}},
+  matrix_of_cc::Vector{BlockSparseMatrix{Matrix{ValType}}},
   rank_of_cc::Vector{Int},
   next_edges_of_cc::Vector{Matrix{Tuple{Vector{Int},Vector{C}}}},
   g::MPOGraph{N,C,Ti},
@@ -162,7 +162,7 @@ On the final site, the scalar `R` factor is applied directly to the local
 blocks.
 """
 @timeit function process_qr(
-  matrix_of_cc::Vector{BlockSparseMatrix{ValType}},
+  matrix_of_cc::Vector{BlockSparseMatrix{Matrix{ValType}}},
   rank_of_cc::Vector{Int},
   next_edges_of_cc::Vector{Matrix{Tuple{Vector{Int},Vector{C}}}},
   g::MPOGraph{N,C,Ti},
diff --git a/src/large-graph-mpo-vertex-cover.jl b/src/large-graph-mpo-vertex-cover.jl
index 67d2dc3..fb56287 100644
--- a/src/large-graph-mpo-vertex-cover.jl
+++ b/src/large-graph-mpo-vertex-cover.jl
@@ -1,12 +1,12 @@
 function _vertex_cover_matrix(
-  ::Type{BlockSparseMatrix{ValType}}, rank::Int
-)::BlockSparseMatrix{ValType} where {ValType}
+  ::Type{BlockSparseMatrix{Matrix{ValType}}}, rank::Int
+)::BlockSparseMatrix{Matrix{ValType}} where {ValType}
   return [Dictionary{Int,Matrix{ValType}}() for _ in 1:rank]
 end
 
 function _vertex_cover_matrix(
-  ::Type{SymbolicBlockSparseMatrix{Ti}}, rank::Int
-)::SymbolicBlockSparseMatrix{Ti} where {Ti<:Integer}
+  ::Type{BlockSparseMatrix{SymbolicLocalMatrix{Ti}}}, rank::Int
+)::BlockSparseMatrix{SymbolicLocalMatrix{Ti}} where {Ti<:Integer}
   return [Dictionary{Int,SymbolicLocalMatrix{Ti}}() for _ in 1:rank]
 end
 
@@ -16,7 +16,7 @@ function _signed_symbolic_local_op_id(lv::LeftVertex, ::Type{Ti})::Ti where {Ti<
 end
 
 function _add_vertex_cover_term!(
-  matrix::BlockSparseMatrix{ValType},
+  matrix::BlockSparseMatrix{Matrix{ValType}},
   m::Int,
   lv::LeftVertex,
   weight::Number,
@@ -32,7 +32,7 @@ function _add_vertex_cover_term!(
 end
 
 function _add_vertex_cover_term!(
-  matrix::SymbolicBlockSparseMatrix{Ti},
+  matrix::BlockSparseMatrix{SymbolicLocalMatrix{Ti}},
   m::Int,
   lv::LeftVertex,
   signed_weight_id::Integer,
@@ -49,7 +49,7 @@ function _add_vertex_cover_term!(
 end
 
 function _add_vertex_cover_unit_term!(
-  matrix::BlockSparseMatrix{ValType},
+  matrix::BlockSparseMatrix{Matrix{ValType}},
   m::Int,
   lv::LeftVertex,
   op_cache::Vector{OpInfo},
@@ -61,7 +61,7 @@ function _add_vertex_cover_unit_term!(
 end
 
 function _add_vertex_cover_unit_term!(
-  matrix::SymbolicBlockSparseMatrix{Ti},
+  matrix::BlockSparseMatrix{SymbolicLocalMatrix{Ti}},
   m::Int,
   lv::LeftVertex,
   op_cache::Vector{OpInfo},
@@ -77,8 +77,8 @@ function _finalize_vertex_cover_matrix!(matrix::BlockSparseMatrix)::Nothing
 end
 
 function _finalize_vertex_cover_matrix!(
-  matrix::SymbolicBlockSparseMatrix
-)::Nothing
+  matrix::BlockSparseMatrix{SymbolicLocalMatrix{Ti}}
+)::Nothing where {Ti}
   for block in matrix
     for terms in values(block)
       _normalize_symbolic_local_matrix!(terms)
diff --git a/src/sparse_tensor_construction.jl b/src/sparse_tensor_construction.jl
index 566ffae..b9241ea 100644
--- a/src/sparse_tensor_construction.jl
+++ b/src/sparse_tensor_construction.jl
@@ -88,7 +88,7 @@ function _fill_splitblocks!(
   blocks::Vector{Block{4}},
   block_values::Vector{C},
   offsets::Vector{Int},
-  block_sparse_matrices::Vector{BlockSparseMatrix{C}},
+  block_sparse_matrices::Vector{BlockSparseMatrix{Matrix{C}}},
 ) where {C}
   entry = 0
 
@@ -125,7 +125,7 @@ are shifted to global right-link ids by `right_link + offsets[i]`.
 """
 function _to_ITensor_splitblocks(
   offsets::Vector{Int},
-  block_sparse_matrices::Vector{BlockSparseMatrix{C}},
+  block_sparse_matrices::Vector{BlockSparseMatrix{Matrix{C}}},
   llink::ITensors.QNIndex,
   rlink::ITensors.QNIndex,
   site::ITensors.QNIndex;
@@ -171,7 +171,7 @@ writes nonzero entries into the appropriate block views.
 """
 function to_ITensor(
   offsets::Vector{Int},
-  block_sparse_matrices::Vector{BlockSparseMatrix{C}},
+  block_sparse_matrices::Vector{BlockSparseMatrix{Matrix{C}}},
   llink::ITensors.QNIndex,
   rlink::ITensors.QNIndex,
   site::ITensors.QNIndex;
@@ -200,11 +200,11 @@ function to_ITensor(
       shifted_right_link = right_link + offset
       right_block = right_blocks[shifted_right_link]
 
-      for (left_link, block) in pairs(column)
+      for (left_link, block_matrix) in pairs(column)
         left_block = left_blocks[left_link]
 
-        @inbounds for site_prime_link in axes(block, 1), site_link in axes(block, 2)
-          value = block[site_prime_link, site_link]
+        @inbounds for site_prime_link in axes(block_matrix, 1), site_link in axes(block_matrix, 2)
+          value = block_matrix[site_prime_link, site_link]
           iszero(value) && continue
 
           push!(
@@ -255,12 +255,12 @@ function to_ITensor(
       right_block = right_blocks[shifted_right_link]
       right_offset = right_offsets[shifted_right_link]
 
-      for (left_link, block) in pairs(column)
+      for (left_link, block_matrix) in pairs(column)
         left_block = left_blocks[left_link]
         left_offset = left_offsets[left_link]
 
-        @inbounds for site_prime_link in axes(block, 1), site_link in axes(block, 2)
-          value = block[site_prime_link, site_link]
+        @inbounds for site_prime_link in axes(block_matrix, 1), site_link in axes(block_matrix, 2)
+          value = block_matrix[site_prime_link, site_link]
           iszero(value) && continue
 
           block_id = _blockid(
@@ -300,7 +300,7 @@ effect for dense indices.
 """
 function to_ITensor(
   offsets::Vector{Int},
-  block_sparse_matrices::Vector{BlockSparseMatrix{C}},
+  block_sparse_matrices::Vector{BlockSparseMatrix{Matrix{C}}},
   llink::Index,
   rlink::Index,
   site::Index;
@@ -311,8 +311,8 @@ function to_ITensor(
     offset = offsets[i]
     matrix = block_sparse_matrices[i]
     for (right_link, column) in enumerate(matrix)
-      for (left_link, block) in pairs(column)
-        tensor[left_link, right_link + offset, :, :] = block
+      for (left_link, block_matrix) in pairs(column)
+        tensor[left_link, right_link + offset, :, :] = block_matrix
       end
     end
   end
diff --git a/test/tensor-construction-tests.jl b/test/tensor-construction-tests.jl
index f9d6d5a..16847bd 100644
--- a/test/tensor-construction-tests.jl
+++ b/test/tensor-construction-tests.jl
@@ -5,7 +5,7 @@ using Test
 
 function reference_to_ITensor(
   offsets::Vector{Int},
-  block_sparse_matrices::Vector{ITensorMPOConstruction.BlockSparseMatrix{C}},
+  block_sparse_matrices::Vector{ITensorMPOConstruction.BlockSparseMatrix{Matrix{C}}},
   llink::ITensors.QNIndex,
   rlink::ITensors.QNIndex,
   site::ITensors.QNIndex;

From 4cd09fe70175add83e9221f06f02c342e400351c Mon Sep 17 00:00:00 2001
From: Ben Corbett <corbett5@llnl.gov>
Date: Fri, 1 May 2026 14:56:21 -0600
Subject: [PATCH 16/16] Added SymbolicWeight

---
 examples/electronic-structure.jl            |  57 +--
 examples/symbolic-transverse-field-ising.jl |   7 +-
 plan.md                                     | 535 --------------------
 src/ITensorMPOConstruction.jl               |   3 +
 src/MPOConstruction.jl                      |  92 +---
 src/OpIDSum.jl                              |  32 +-
 src/SymbolicLocalMatrix.jl                  | 154 +-----
 src/SymbolicWeight.jl                       |  28 +
 src/large-graph-mpo-common.jl               |  19 +-
 src/large-graph-mpo-vertex-cover.jl         |  68 +--
 test/symbolic-local-matrix-tests.jl         | 185 ++-----
 test/symbolic-mpo-tests.jl                  |  57 +--
 12 files changed, 181 insertions(+), 1056 deletions(-)
 delete mode 100644 plan.md
 create mode 100644 src/SymbolicWeight.jl

diff --git a/examples/electronic-structure.jl b/examples/electronic-structure.jl
index 72cba6b..6e7c1d9 100644
--- a/examples/electronic-structure.jl
+++ b/examples/electronic-structure.jl
@@ -206,16 +206,6 @@ end
 
 # # Symbolic construction
 
-function electronic_structure_fermion_sort_sign(site_ids::Integer...)::Int
-  sign = 1
-  for i in eachindex(site_ids)
-    for j in 1:(i - 1)
-      site_ids[i] < site_ids[j] && (sign = -sign)
-    end
-  end
-  return sign
-end
-
 function electronic_structure_symbolic_OpIDSum(
   N::Int,
   op_cache_vec::OpCacheVec,
@@ -226,16 +216,15 @@ function electronic_structure_symbolic_OpIDSum(
   opC(k::Int, spin::Bool) = OpID(2 + spin, k)
   opCdag(k::Int, spin::Bool) = OpID(4 + spin, k)
 
-  os = OpIDSum{4,Int,Int}(2 * N^4 + 2 * N^2, op_cache_vec)
+  os = OpIDSum{4,SimpleWeight,Int}(2 * N^4 + 2 * N^2, op_cache_vec)
 
   map_1e = Tuple{Int, Int}[]
   for p in 1:N
     for q in 1:N
       push!(map_1e, (p, q))
       id = length(map_1e)
-      sign = electronic_structure_fermion_sort_sign(p, q)
       for spin in (↓, ↑)
-        add!(os, sign * id, opCdag(p, spin), opC(q, spin))
+        add!(os, SimpleWeight(id), opCdag(p, spin), opC(q, spin))
       end
     end
   end
@@ -260,10 +249,9 @@ function electronic_structure_symbolic_OpIDSum(
             push!(map_2e, (p, q, r, s, s_p, s_q, s_r, s_s))
 
             id = length(map_2e) + length(map_1e)
-            sign = electronic_structure_fermion_sort_sign(p, q, r, s)
             add!(
               os,
-              sign * id,
+              SimpleWeight(id),
               opCdag(p, s_p),
               opCdag(q, s_q),
               opC(r, s_r),
@@ -289,7 +277,7 @@ function electronic_structure_symbolic_coefficients(
   )
 
   for (id, (p, q)) in pairs(map_1e)
-    coefficients[id] = electronic_structure_fermion_sort_sign(p, q) * h[p, q]
+    coefficients[id] = h[p, q]
   end
 
   offset = length(map_1e)
@@ -308,7 +296,7 @@ function electronic_structure_symbolic_coefficients(
       coeff += V[q, r, p, s]
     end
 
-    coefficients[offset + i] = electronic_structure_fermion_sort_sign(p, q, r, s) * coeff
+    coefficients[offset + i] = coeff
   end
 
   return coefficients
@@ -320,7 +308,7 @@ function mpo_relative_difference(A::MPO, B::MPO)::Float64
 end
 
 let
-  N = 50
+  N = 6
   println("Constructing a symbolic electronic structure MPO for $N sites")
 
   h, V = get_coefficients(N)
@@ -356,26 +344,19 @@ let
 
   print_timer()
 
-  # @time "Constructing reference MPO" H_numeric = MPO_new(
-  #   numeric_os,
-  #   sites;
-  #   alg="VC",
-  #   basis_op_cache_vec=op_cache_vec,
-  #   splitblocks=true,
-  #   check_for_errors=false,
-  #   checkflux=false,
-  # )
-
-  # relative_difference = mpo_relative_difference(H_symbolic, H_numeric)
-  # println("Symbolic instantiation relative difference: $relative_difference")
-  # @assert relative_difference < 1e-10
-
-  # h2 = 1.1 .* h
-  # V2 = 1.1 .* V
-  # coefficients2 = electronic_structure_symbolic_coefficients(h2, V2, map_1e, map_2e)
-  # H_symbolic2 = instantiate_MPO(sym, coefficients2; splitblocks=true, checkflux=false)
-  # @assert mpo_relative_difference(H_symbolic2, H_symbolic) > 0
-  # println("Updated symbolic coefficients without rebuilding the symbolic MPO")
+  @time "Constructing reference MPO" H_numeric = MPO_new(
+    numeric_os,
+    sites;
+    alg="VC",
+    basis_op_cache_vec=op_cache_vec,
+    splitblocks=true,
+    check_for_errors=false,
+    checkflux=false,
+  )
+
+  relative_difference = mpo_relative_difference(H_symbolic, H_numeric)
+  println("Symbolic instantiation relative difference: $relative_difference")
+  @assert relative_difference < 1e-10
 end
 
 nothing;
diff --git a/examples/symbolic-transverse-field-ising.jl b/examples/symbolic-transverse-field-ising.jl
index 9168b33..03e92ed 100644
--- a/examples/symbolic-transverse-field-ising.jl
+++ b/examples/symbolic-transverse-field-ising.jl
@@ -24,19 +24,18 @@ function all_to_all_tfim_symbolic(N::Int)::SymbolicMPO
   Z(n::Int) = OpID(3, n)
 
   num_couplings = num_tfim_couplings(N)
-  field_label = num_couplings + 1
-  os = OpIDSum{2,Int,Int}(num_couplings + N, op_cache_vec)
+  os = OpIDSum{2,SimpleWeight,Int}(num_couplings + N, op_cache_vec)
 
   label = 1
   for i in 1:N
     for j in (i + 1):N
-      add!(os, label, Z(i), Z(j))
+      add!(os, SimpleWeight(label), Z(i), Z(j))
       label += 1
     end
   end
 
   for i in 1:N
-    add!(os, field_label, X(i))
+    add!(os, SimpleWeight(label), X(i))
   end
 
   return MPO_symbolic(os, sites)
diff --git a/plan.md b/plan.md
deleted file mode 100644
index 59d6e53..0000000
--- a/plan.md
+++ /dev/null
@@ -1,535 +0,0 @@
-# Symbolic Vertex-Cover MPO Construction
-
-## Summary
-
-Add a reusable symbolic MPO construction path for integer-labeled `OpIDSum`s.
-`MPO_symbolic` should run the vertex-cover MPO construction once, store the
-resulting tensor entries as integer linear combinations of coefficient
-variables, and allow repeated numeric instantiation with different coefficient
-vectors.
-
-The primary API should follow the package's current argument order:
-
-```julia
-sym = MPO_symbolic(os, sites; basis_op_cache_vec=nothing,
-                   check_for_errors=true,
-                   combine_qn_sectors=false,
-                   output_level=0)
-
-H = instantiate_MPO(sym, coefficients; splitblocks=true, checkflux=false)
-instantiate_MPO!(H, sym, coefficients; checkflux=false)
-H2 = instantiate_MPO(H, sym, coefficients; checkflux=false)
-```
-
-## Public API and V1 Semantics
-
-- Export `SymbolicMPO` and `MPO_symbolic`.
-- Support only `OpIDSum{N,<:Integer,Ti}` input in V1.
-- Symbolic construction is vertex-cover construction. It should not expose an
-  `alg` keyword, and there is no QR symbolic path to support.
-- Do not add symbolic `OpSum` conversion in V1.
-- Preserve numeric `MPO_new` behavior exactly.
-- User-provided integer labels are variable labels and must always be greater
-  than zero.
-- Convert each user label `k` to internal symbolic id `k + 1`, reserving
-  internal id `1` for the constant one introduced by the vertex-cover algorithm.
-- Carry signs by negating the internal symbolic id. Negative ids can arise from
-  sign-only basis rewrites and other internal sign handling; they are not
-  user-facing labels.
-- During numeric substitution:
-  - internal id `+1` evaluates to `+1.0`,
-  - internal id `-1` evaluates to `-1.0`,
-  - internal id `k > 1` evaluates to `coefficients[k - 1]`,
-  - internal id `k < -1` evaluates to `-coefficients[-k - 1]`.
-- The coefficient vector is indexed by the original positive user labels. The
-  reserved internal constant id is transparent to the user, so there is no
-  reserved slot in `coefficients`.
-- Reject `OpSum` input, unsupported basis rewrite factors, and missing
-  coefficients with clear `ArgumentError`s.
-
-## Internal Representation
-
-### Symbolic local matrix terms
-
-Do not model symbolic MPO entries as scalar coefficients. A symbolic block entry
-is a symbolic weighted sum of local operator matrices:
-
-```julia
-const SymbolicLocalMatrix{Ti<:Integer} = Vector{Tuple{Int,Ti}}
-```
-
-Each tuple is `(signed_weight_id, signed_local_op_id)`.
-
-- `signed_weight_id` is the signed internal coefficient id:
-  - `+1` and `-1` are the constants `+1` and `-1`,
-  - `k > 1` maps to user assignment `coefficients[k - 1]`,
-  - `k < -1` maps to `-coefficients[-k - 1]`.
-- `abs(signed_local_op_id)::Ti` uses the same integer type as
-  `OpIDSum{N,C,Ti}` and identifies the local operator being scaled by that
-  signed coefficient id.
-- `signed_local_op_id < 0` means the local operator should be emitted with the
-  Jordan-Wigner string applied. `signed_local_op_id > 0` means no
-  Jordan-Wigner string is applied.
-- A term represents
-  `substitute(signed_weight_id, coefficients) * local_matrix`, where
-  `local_matrix` is `op_cache_vec[n][abs(signed_local_op_id)].matrix`, with the
-  Jordan-Wigner sign pattern applied when `signed_local_op_id < 0`.
-
-Required helpers:
-
-- `_weight_value(signed_weight_id, coefficients)` substitutes one signed
-  internal id using the rules above.
-- `_normalize_symbolic_local_matrix!(terms)` canonicalizes ordering and removes
-  canceling pairs such as `(k, op)` with `(-k, op)`. It must not collapse
-  repeated identical pairs into one entry, since multiplicity represents
-  repeated contributions such as `2 * coefficients[k - 1] * op`.
-- `_scale_symbolic_weight(signed_weight_id, factor)` accepts only exact factors
-  `+1` and `-1` in V1, flipping the signed id for `-1`.
-- `_evaluate_symbolic_local_matrix(terms, coefficients, op_cache)` builds the
-  dense numeric local matrix by summing cached local operator matrices, applying
-  Jordan-Wigner strings for negative local operator ids, and substituting
-  weights.
-
-### Symbolic Block Storage
-
-Do not reuse `BlockSparseMatrix{C}` for symbolic block entries, since that alias
-stores dense `Matrix{C}` blocks. Symbolic construction should use:
-
-```julia
-const SymbolicBlockSparseMatrix{Ti<:Integer} =
-  Vector{Dictionary{Int,SymbolicLocalMatrix{Ti}}}
-```
-
-This mirrors the existing `right_link => left_link => block` layout, but each
-block is a symbolic local matrix rather than an already-expanded dense matrix.
-Numeric instantiation converts these symbolic blocks into ordinary
-`BlockSparseMatrix{T}` values immediately before constructing `ITensor`s.
-
-## Implementation Checklist
-
-### 1. Add and test symbolic local matrix operations
-
-- Add `SymbolicLocalMatrix{Ti}` and `SymbolicBlockSparseMatrix{Ti}` aliases or
-  small wrapper structs.
-- Implement helpers for appending terms, normalizing terms, applying `+1` or
-  `-1` rewrite factors, computing `max_user_label`, and substituting signed
-  internal ids.
-- Implement `_evaluate_symbolic_local_matrix` to produce a dense numeric matrix
-  for one site from `(signed_weight_id, signed_local_op_id)` terms and the site's
-  operator cache.
-- Add focused tests for:
-  - duplicate `(weight_id, signed_local_op_id)` terms preserving multiplicity,
-  - zero normalization,
-  - `+1` and `-1` constant substitution,
-  - positive user label to internal id mapping,
-  - signed internal id evaluation,
-  - missing coefficient rejection,
-  - non-`±1` rewrite-factor rejection.
-
-Tests for this step:
-
-- Unit-test `_weight_value` directly:
-  - `+1` returns `1`,
-  - `-1` returns `-1`,
-  - `+3` returns `coefficients[2]`,
-  - `-3` returns `-coefficients[2]`,
-  - `+4` throws when `length(coefficients) == 2`.
-- Unit-test symbolic local matrix normalization:
-  - `[(3, 2), (-3, 2)]` normalizes to an empty term list,
-  - `[(3, 2), (3, 2)]` keeps two entries and evaluates to twice the same local
-    matrix contribution,
-  - terms are sorted deterministically so equality checks are stable.
-- Unit-test Jordan-Wigner encoding:
-  - `(3, -2)` evaluates to `coefficients[2]` times operator `2` with the
-    Jordan-Wigner sign pattern,
-  - `(3, 2)` evaluates to the same operator without the Jordan-Wigner sign
-    pattern.
-
-### 2. Convert the integer-labeled `OpIDSum` into internal-id form
-
-- Add an internal conversion helper:
-
-```julia
-internalize_symbolic_ids!(os::OpIDSum{N,C,Ti}) where {N,C<:Integer,Ti}
-```
-
-- This helper may mutate the provided `OpIDSum`, matching the existing MPO
-  construction behavior where preprocessing, sorting, and duplicate compaction
-  are destructive.
-- Map each stored positive scalar label to an internal id in place.
-- Reject scalar labels less than or equal to `0`.
-- Map stored scalar label `k > 0` to internal id `k + 1`.
-- Keep the same `op_cache_vec`, `abs_tol`, and `modify!`.
-- Add tests proving the in-place conversion maps labels correctly and rejects
-  nonpositive labels.
-
-Tests for this step:
-
-- Build a small `OpIDSum{2,Int,Int}` with scalar labels `1`, `2`, and `3`;
-  after `internalize_symbolic_ids!`, assert that the stored labels are `2`, `3`,
-  and `4`.
-- Add terms with labels `0` and `-1` and assert `internalize_symbolic_ids!`
-  throws an `ArgumentError`.
-- Assert the helper is intentionally mutating by checking that the same `os`
-  object has updated scalar storage after conversion.
-- Assert `op_cache_vec`, `abs_tol`, and `modify!` are still the original objects
-  or equivalent values after conversion.
-
-### 3. Make basis rewrite compatible with internal signed ids
-
-- Reuse `prepare_opID_sum!(os, basis_op_cache_vec; symbolic_coefficients=true)`
-  after the in-place internal-id conversion. The keyword keeps ordinary numeric
-  `MPO_new` preprocessing on the existing coefficient-scaling path while the
-  symbolic path treats integer coefficients as internal signed ids.
-- Ensure `rewrite_in_operator_basis!` can multiply internal signed ids by basis
-  rewrite factors only when the factor is exactly `+1` or `-1`.
-- Reject `im`, `2im`, `0.5`, `2`, and other unsupported factors with an
-  `ArgumentError`; this may be reported through Julia's threaded task failure
-  wrapper because symbolic rewrites use the same threaded loop as numeric
-  rewrites.
-- Add a sign-only rewrite regression where a same-site product produces `-1`
-  and verifies the internal id sign flips while the user label is preserved.
-
-Tests for this step:
-
-- Construct a small operator-cache fixture where a same-site product rewrites
-  into a basis operator with factor `-1`; after preparation, verify internal id
-  `k + 1` becomes `-(k + 1)` and the basis operator id is updated.
-- Construct a fixture with rewrite factor `+1`; verify the internal id sign is
-  unchanged.
-- Reuse the existing `X * Y -> im * Z` style case or an equivalent fixture to
-  assert complex factors throw.
-- Add separate fixtures for `0.5` and `2` factors and assert both throw, since
-  V1 accepts only exact `±1`.
-
-### 4. Build a symbolic graph without adding identifier integers
-
-- Numeric `MPOGraph(os::OpIDSum)` must keep its current duplicate-term
-  compaction and tolerance behavior.
-- Add a symbolic graph-construction path for the internal-id `OpIDSum`.
-- Duplicate operator tuples must combine correctly even when their scalar labels
-  are different symbolic variables. Combining duplicates should preserve all
-  signed ids, for example by storing a vector of signed internal ids as the edge
-  weight or by expanding the duplicates into distinct edge entries before local
-  block construction.
-- The graph edge weight type should remain a signed internal id or a collection
-  of signed ids, not a symbolic local matrix term.
-- Add a test where duplicate operator tuples with different user labels both
-  contribute to the same symbolic block and instantiate correctly.
-
-Tests for this step:
-
-- Build an `OpIDSum{2,Int,Int}` with two identical operator tuples labeled `1`
-  and `2`; after symbolic graph construction and duplicate-right-vertex
-  compaction, assert the resulting graph preserves both signed ids as separate
-  edge weights, not their numeric sum.
-- Add two identical operator tuples with the same label `1`; instantiate with
-  `coefficients[1] = a` through the symbolic local matrix helper and verify the
-  local contribution is `2a`.
-- Add two identical operator tuples whose labels become opposite signed internal
-  ids after preprocessing; verify the derived symbolic local matrix cancels or
-  evaluates to zero for that contribution.
-
-### 5. Build `SymbolicMPO`
-
-Add:
-
-```julia
-struct SymbolicMPO{Ti<:Integer}
-  offsets::Vector{Vector{Int}}
-  block_sparse_matrices::Vector{Vector{SymbolicBlockSparseMatrix{Ti}}}
-  sites::Vector{<:Index}
-  llinks::Vector{Index}
-  op_cache_vec::OpCacheVec
-  max_user_label::Int
-end
-```
-
-If Julia requires concrete field types for `sites`, parameterize the site vector
-type instead of using an abstract field.
-
-`MPO_symbolic(os, sites; kwargs...)` should:
-
-- validate `os` has an integer coefficient type,
-- reject any supplied `alg` keyword with an `ArgumentError` explaining that
-  symbolic construction is always vertex-cover based,
-- convert `basis_op_cache_vec` with `to_OpCacheVec(sites, basis_op_cache_vec)`,
-- convert `os` to internal-id form in place,
-- run `prepare_opID_sum!` and optionally `check_os_for_errors`,
-- build the symbolic graph without numerically adding identifier ids,
-- initialize the dummy left link exactly as `MPO_new` does,
-- call the existing `resume_MPO_construction!` with `alg="VC"` and symbolic
-  block storage, so symbolic construction reuses `at_site!` and
-  `process_vertex_cover!` instead of maintaining parallel symbolic driver
-  functions,
-- use small storage-dispatch helpers inside the vertex-cover path so numeric
-  blocks still accumulate dense matrices while symbolic blocks append
-  `(signed_weight_id, signed_local_op_id)` terms,
-- when adding a local contribution, append
-  `(signed_weight_id, signed_local_op_id)` to the symbolic block instead of
-  expanding the dense local operator matrix,
-- store offsets, symbolic block matrices, sites, links, `op_cache_vec`, and the
-  maximum positive user label needed for substitution.
-
-Tests for this step:
-
-- Construct a small transverse-field Ising `OpIDSum{2,Int,Int}` and assert
-  `MPO_symbolic(os, sites)` returns a `SymbolicMPO` with:
-  - `length(offsets) == length(sites)`,
-  - `length(block_sparse_matrices) == length(sites)`,
-  - `length(llinks) == length(sites) + 1`,
-  - `max_user_label` equal to the largest user label in the input.
-- Assert `MPO_symbolic(os, sites; alg="QR")` throws an `ArgumentError`.
-- Assert `MPO_symbolic(::OpSum, sites)` throws an `ArgumentError` or has no
-  method, whichever API choice is implemented.
-- For a fermionic fixture, inspect at least one symbolic local matrix term and
-  assert negative `signed_local_op_id` is used where `needs_JW_string` is true.
-
-### 6. Instantiate a fresh numeric MPO
-
-Add:
-
-```julia
-instantiate_MPO(sym::SymbolicMPO, coefficients; splitblocks=true, checkflux=false)
-```
-
-It should:
-
-- validate the coefficient vector length against `sym.max_user_label`,
-- at each site, evaluate `sym.block_sparse_matrices[n]` with `sym.op_cache_vec`
-  into a `Vector{BlockSparseMatrix{C}}`,
-- preserve the symbolic link layout even when evaluated values are zero,
-- call the existing tensor assembly path with the requested `splitblocks` value,
-- contract away the dummy boundary links the same way as the low-level
-  `instantiate_MPO(offsets, block_sparse_matrices, sites, llinks; ...)`,
-- optionally run `ITensors.checkflux`.
-
-For QN tensors, structural block discovery uses the existing numeric
-`to_ITensor` path after symbolic block evaluation. A coefficient vector
-containing zeros must not change the stored symbolic MPO link layout.
-
-Tests for this step:
-
-- Instantiate a small symbolic MPO with two different coefficient vectors and
-  compare each result against numeric `MPO_new` built from an equivalent numeric
-  `OpIDSum`.
-- Run the duplicate-label symbolic graph regression through `MPO_symbolic` and
-  compare the instantiated MPO against numeric `MPO_new` built from the
-  corresponding numeric coefficients.
-- Instantiate with `splitblocks=false` and `splitblocks=true` for a QN-conserving
-  fixture and compare both against numeric `MPO_new`.
-- Instantiate with a coefficient vector containing zeros and assert the link
-  dimensions match a nonzero-coefficient instantiation from the same
-  `SymbolicMPO`.
-- Pass a too-short coefficient vector and assert an `ArgumentError`.
-- Run with `checkflux=true` on a QN fixture and assert it succeeds.
-
-### 7. Add in-place and template-assisted instantiation
-
-Add:
-
-```julia
-instantiate_MPO!(H::MPO, sym::SymbolicMPO, coefficients; checkflux=false)
-instantiate_MPO(H_template::MPO, sym::SymbolicMPO, coefficients; checkflux=false)
-```
-
-These overloads must not accept a `splitblocks` keyword. A preconstructed MPO
-already fixes the block layout, so evaluation must use the layout of `H` or
-`H_template`.
-
-`instantiate_MPO!` should:
-
-- check `length(H) == length(sym.sites)`,
-- check site indices are compatible with `sym.sites`,
-- check link dimensions match `sym.llinks`,
-- temporarily add the dummy boundary links back to the first and last MPO
-  tensors by multiplying by one-entry dummy tensors,
-- for each site, zero the existing entries with `H[n] .= 0` and refill through
-  the same four-index coordinate layout used by the numeric tensor assembly,
-- contract away the dummy boundary links again at the end,
-- preserve the existing `MPO` object identity and link layout.
-
-`instantiate_MPO(H_template, sym, coefficients)` should:
-
-- reject incompatible templates with `ArgumentError`,
-- copy a compatible template's tensors before filling them,
-- call `instantiate_MPO!`,
-- return the instantiated MPO.
-
-This path should avoid rebuilding ITensors for compatible templates; the
-expected fast path is to reuse the existing structural storage for the nonzero
-entries after temporarily restoring the boundary-link layout.
-
-Tests for this step:
-
-- Build `H = instantiate_MPO(sym, coeffs1; splitblocks=true)`, then call
-  `instantiate_MPO!(H, sym, coeffs2)` and compare `H` against
-  `instantiate_MPO(sym, coeffs2; splitblocks=true)`.
-- Assert the `MPO` object identity and link dimensions are unchanged after
-  `instantiate_MPO!`.
-- Assert `instantiate_MPO!(H, sym, coeffs; splitblocks=false)` throws a
-  `MethodError` or `ArgumentError`, since this overload must not accept
-  `splitblocks`.
-- Call `instantiate_MPO(H, sym, coeffs2)` with a compatible template and compare
-  against fresh instantiation using the same layout.
-- Assert template-assisted instantiation does not mutate the template MPO.
-- Build an incompatible template with different sites or link dimensions and
-  assert `instantiate_MPO(H_bad, sym, coeffs)` throws an `ArgumentError`.
-
-### 8. Add symbolic MPO tests
-
-Create `test/symbolic-mpo-tests.jl` and include it from `test/runtests.jl`.
-
-The file should use named `@testset`s matching the implementation steps, so a
-failure points directly to the layer that broke:
-
-- `"symbolic local matrix terms"`:
-  - normalization,
-  - duplicate `(weight_id, signed_local_op_id)` multiplicity,
-  - constant weight ids,
-  - negative local operator id applies the Jordan-Wigner string,
-  - positive user label to internal id mapping,
-  - signed internal ids,
-  - missing coefficient error,
-  - unsupported rewrite-factor error.
-- `"internal symbolic ids"`:
-  - in-place user-label to internal-id conversion,
-  - nonpositive label rejection,
-  - negative internal ids from later sign-carrying operations remain supported.
-- `"basis rewrite"`:
-  - `+1` rewrite factor preserves id sign,
-  - `-1` rewrite factor flips id sign,
-  - complex, fractional, and integer factors other than `±1` throw.
-- `"symbolic graph duplicate terms"`:
-  - different labels on the same operator tuple are both preserved,
-  - repeated identical labels preserve multiplicity,
-  - opposite signed ids cancel or instantiate to zero.
-- `"fresh instantiation"`:
-  - build `OpIDSum{2,Int,Int}`,
-  - construct `sym = MPO_symbolic(os, sites)`,
-  - instantiate with several random coefficient vectors indexed directly by the
-    positive user labels, passing `splitblocks` to fresh `instantiate_MPO`,
-  - compare against numeric `MPO_new` built with the same coefficients.
-- `"public API rejections"`:
-  - any supplied `alg` keyword,
-  - `OpSum` input,
-  - missing coefficient,
-  - incompatible preconstructed MPO templates.
-- `"in-place instantiation"`:
-  - instantiate once,
-  - update with a second coefficient vector through `instantiate_MPO!`,
-  - verify `instantiate_MPO!` does not accept a `splitblocks` keyword,
-  - compare to fresh instantiation,
-  - check link dimensions remain unchanged.
-- `"template-assisted instantiation"`:
-  - instantiate from a compatible template,
-  - compare to fresh instantiation,
-  - reject an incompatible template.
-
-Tests for this step:
-
-- `julia --project test/symbolic-mpo-tests.jl` runs all new symbolic testsets
-  by itself.
-- `julia --project test/runtests.jl` includes `test/symbolic-mpo-tests.jl` once
-  and still runs the existing test files.
-- Existing `test/ops-tests.jl` still passes, proving the basis-rewrite changes
-  did not regress the numeric `OpIDSum` path.
-
-### 9. Update electronic-structure example
-
-Add a compact symbolic section to `examples/electronic-structure.jl`.
-
-Implementation outline:
-
-- Build a small integer-labeled `OpIDSum` in the same term order as the numeric
-  electronic-structure `OpIDSum`.
-- Keep a `coefficients` vector indexed directly by positive user labels.
-- Whenever user label `k` is added to the symbolic `OpIDSum`, store the actual
-  numeric coefficient at `coefficients[k]`.
-- Construct `sym = MPO_symbolic(
-  os_symbolic, sites; basis_op_cache_vec=os_symbolic.op_cache_vec, check_for_errors=false
-  )`.
-- Instantiate twice with two coefficient vectors, passing `splitblocks=true` to
-  fresh `instantiate_MPO`.
-- For a small `N`, compare one symbolic instantiation to the existing numeric
-  `MPO_new` result.
-- Keep the example section short so the existing benchmark narrative remains
-  readable.
-
-Tests for this step:
-
-- Run the example with a small `N` setting first and assert the symbolic
-  instantiation matches the numeric `MPO_new` result using the existing MPO
-  comparison helper or an equivalent relative-difference check.
-- Run the final example command:
-
-```bash
-julia --project -t8 examples/electronic-structure.jl
-```
-
-- Confirm the example still reports the existing numeric construction metrics
-  and additionally prints or checks the symbolic construction section.
-- Confirm changing the second coefficient vector changes the instantiated MPO
-  without rebuilding the symbolic structure.
-
-### 10. Update docs
-
-- Add docstrings for `SymbolicMPO`, `MPO_symbolic`, and the new
-  `instantiate_MPO` / `instantiate_MPO!` overloads.
-- Update `docs/src/documentation/MPO_new.md` to include the symbolic API in the
-  `@docs` block.
-- Add a short prose section explaining:
-  - integer labels in `OpIDSum`,
-  - positive user labels and transparent internal remapping,
-  - repeated instantiation,
-  - why symbolic construction is vertex-cover based and has no QR option,
-  - `±1` rewrite-factor limitation,
-  - `OpIDSum`-only limitation.
-
-Tests for this step:
-
-- Verify `julia --project -e 'using ITensorMPOConstruction'` loads cleanly after
-  adding exports and docstrings.
-- Verify the docs build:
-
-```bash
-julia --project=docs --color=yes docs/make.jl
-```
-
-- Check the generated docs page includes `SymbolicMPO`, `MPO_symbolic`,
-  `instantiate_MPO(sym, ...)`, and `instantiate_MPO!(H, sym, ...)`.
-
-## Verification Checklist
-
-Run focused checks first:
-
-```bash
-julia --project -e 'using ITensorMPOConstruction'
-julia --project test/symbolic-mpo-tests.jl
-julia --project test/ops-tests.jl
-julia --project test/test-MPOConstruction.jl
-```
-
-Then run broader checks:
-
-```bash
-julia --project test/runtests.jl
-julia --project=docs --color=yes docs/make.jl
-```
-
-Because `examples/electronic-structure.jl` is modified, also run:
-
-```bash
-julia --project -t8 examples/electronic-structure.jl
-```
-
-## Acceptance Criteria
-
-- Existing numeric `MPO_new` tests pass unchanged.
-- `MPO_symbolic` may mutate its input `OpIDSum`, consistently with existing MPO
-  construction preprocessing.
-- Repeated symbolic instantiation matches fresh numeric MPO construction.
-- In-place symbolic instantiation updates an existing MPO without changing link
-  dimensions.
-- Symbolic construction rejects unsupported input cases with clear errors.
-- Docs build successfully and include the new public API.
diff --git a/src/ITensorMPOConstruction.jl b/src/ITensorMPOConstruction.jl
index 836fa4b..d73ffb9 100644
--- a/src/ITensorMPOConstruction.jl
+++ b/src/ITensorMPOConstruction.jl
@@ -18,6 +18,7 @@ using ThreadsX
 #
 include("time_if.jl")
 include("OpIDSum.jl")
+include("SymbolicWeight.jl")
 include("SymbolicLocalMatrix.jl")
 include("ops.jl")
 include("BipartiteGraph.jl")
@@ -36,6 +37,8 @@ include("MPOConstruction.jl")
 # OpIDSum.jl
 export OpInfo, OpCacheVec, to_OpCacheVec, OpID, OpIDSum
 
+export SimpleWeight
+
 # MPOConstruction.jl
 export resume_MPO_construction!,
   MPO_new, SymbolicMPO, MPO_symbolic, instantiate_MPO, instantiate_MPO!, sparsity, block2_nnz
diff --git a/src/MPOConstruction.jl b/src/MPOConstruction.jl
index cff29e8..b4523c0 100644
--- a/src/MPOConstruction.jl
+++ b/src/MPOConstruction.jl
@@ -45,7 +45,7 @@ $_resume_MPO_kwargs
 function resume_MPO_construction!(
   n_init::Int,
   offsets::Vector{Vector{Int}},
-  block_sparse_matrices::Vector{Vector{MatrixType}},
+  block_sparse_matrices::Vector{Vector{BlockSparseMatrix{MatrixType}}},
   sites::Vector{<:Index},
   llinks::Vector{<:Index},
   g::MPOGraph,
@@ -57,7 +57,7 @@ function resume_MPO_construction!(
   call_back::Function=(
     cur_site::Int,
     offsets::Vector{Vector{Int}},
-    block_sparse_matrices::Vector{Vector{MatrixType}},
+    block_sparse_matrices::Vector{Vector{BlockSparseMatrix{MatrixType}}},
     sites::Vector{<:Index},
     llinks::Vector{<:Index},
     cur_graph::MPOGraph,
@@ -159,24 +159,13 @@ MPO construction, but each block entry is a symbolic local matrix term list.
 later numeric substitution.
 """
 struct SymbolicMPO{
-  Ti<:Integer,Sites<:AbstractVector{<:Index},Links<:AbstractVector{<:Index}
+  W<:AbstractWeight,Ti<:Integer,Sites<:AbstractVector{<:Index},Links<:AbstractVector{<:Index}
 }
   offsets::Vector{Vector{Int}}
-  block_sparse_matrices::Vector{Vector{BlockSparseMatrix{SymbolicLocalMatrix{Ti}}}}
+  block_sparse_matrices::Vector{Vector{BlockSparseMatrix{SymbolicLocalMatrix{W, Ti}}}}
   sites::Sites
   llinks::Links
   op_cache_vec::OpCacheVec
-  max_user_label::Int
-end
-
-function _check_symbolic_coefficients(sym::SymbolicMPO, coefficients::AbstractVector)::Nothing
-  length(coefficients) == sym.max_user_label || throw(
-    ArgumentError(
-      "Expected $(sym.max_user_label) coefficients. " *
-      "Received $(length(coefficients)) coefficients.",
-    ),
-  )
-  return nothing
 end
 
 function _symbolic_mpo_eltype(sym::SymbolicMPO, coefficients::AbstractVector)::Type
@@ -190,16 +179,16 @@ function _symbolic_mpo_eltype(sym::SymbolicMPO, coefficients::AbstractVector)::T
 end
 
 function _evaluate_symbolic_block_sparse_matrix(
-  symbolic_matrix::BlockSparseMatrix{SymbolicLocalMatrix{Ti}},
+  symbolic_matrix::BlockSparseMatrix{SymbolicLocalMatrix{W, Ti}},
   coefficients::AbstractVector,
   op_cache::Vector{OpInfo},
-  ::Type{C},
-)::BlockSparseMatrix{Matrix{C}} where {Ti,C}
-  matrix = [Dictionary{Int,Matrix{C}}() for _ in eachindex(symbolic_matrix)]
+  ::Type{ValType},
+)::BlockSparseMatrix{Matrix{ValType}} where {W,Ti,ValType}
+  matrix = BlockSparseMatrix(Matrix{ValType}, length(symbolic_matrix))
 
   for right_link in eachindex(symbolic_matrix)
     for (left_link, terms) in pairs(symbolic_matrix[right_link])
-      block = _evaluate_symbolic_local_matrix(C, terms, coefficients, op_cache)
+      block = _evaluate_symbolic_local_matrix(ValType, terms, coefficients, op_cache)
       set!(matrix[right_link], left_link, block)
     end
   end
@@ -208,15 +197,15 @@ function _evaluate_symbolic_block_sparse_matrix(
 end
 
 function _evaluate_symbolic_block_sparse_matrices(
-  symbolic_matrices::Vector{BlockSparseMatrix{SymbolicLocalMatrix{Ti}}},
+  symbolic_matrices::Vector{BlockSparseMatrix{SymbolicLocalMatrix{W, Ti}}},
   coefficients::AbstractVector,
   op_cache::Vector{OpInfo},
-  ::Type{C},
-)::Vector{BlockSparseMatrix{Matrix{C}}} where {Ti,C}
-  matrices = Vector{BlockSparseMatrix{Matrix{C}}}(undef, length(symbolic_matrices))
+  ::Type{ValType},
+)::Vector{BlockSparseMatrix{Matrix{ValType}}} where {W, Ti,ValType}
+  matrices = Vector{BlockSparseMatrix{Matrix{ValType}}}(undef, length(symbolic_matrices))
   for cc in eachindex(symbolic_matrices)
     matrices[cc] = _evaluate_symbolic_block_sparse_matrix(
-      symbolic_matrices[cc], coefficients, op_cache, C
+      symbolic_matrices[cc], coefficients, op_cache, ValType
     )
   end
   return matrices
@@ -236,7 +225,6 @@ assembly path used by numeric MPO construction.
 function instantiate_MPO(
   sym::SymbolicMPO, coefficients::AbstractVector; splitblocks::Bool=true, checkflux::Bool=false
 )::MPO
-  _check_symbolic_coefficients(sym, coefficients)
   C = _symbolic_mpo_eltype(sym, coefficients)
 
   H = MPO(sym.sites)
@@ -257,15 +245,7 @@ function instantiate_MPO(
     )
   end
 
-  # Remove the dummy link indices from the left and right.
-  L = ITensor(sym.llinks[1])
-  L[end] = 1.0
-
-  R = ITensor(dag(sym.llinks[end]))
-  R[1] = 1.0
-
-  H[1] *= L
-  H[length(sym.sites)] *= R
+  _remove_symbolic_mpo_boundary_links!(H, sym.llinks)
 
   if checkflux
     @timeit "checkflux" Threads.@threads for n in eachindex(sym.sites)
@@ -551,7 +531,6 @@ object `H` is updated in place and returned.
 function instantiate_MPO!(
   H::MPO, sym::SymbolicMPO, coefficients::AbstractVector; checkflux::Bool=false
 )::MPO
-  _check_symbolic_coefficients(sym, coefficients)
   _check_symbolic_mpo_template(H, sym)
 
   _add_symbolic_mpo_boundary_links!(H, sym)
@@ -596,22 +575,6 @@ function instantiate_MPO(
   return instantiate_MPO!(deepcopy(H_template), sym, coefficients; checkflux)
 end
 
-function _check_MPO_symbolic_kwargs(kwargs)::Nothing
-  if haskey(kwargs, :alg)
-    throw(
-      ArgumentError(
-        "Symbolic MPO construction is always vertex-cover based and does not accept an alg keyword.",
-      ),
-    )
-  end
-
-  if !isempty(kwargs)
-    throw(ArgumentError("Unsupported keyword(s) for MPO_symbolic: $(join(keys(kwargs), ", "))."))
-  end
-
-  return nothing
-end
-
 """
     MPO_symbolic(os::OpIDSum, sites; basis_op_cache_vec=nothing,
       check_for_errors=true, combine_qn_sectors=false, output_level=0) -> SymbolicMPO
@@ -624,27 +587,19 @@ coefficient labels are remapped to internal symbolic ids, an optional basis
 rewrite is applied, and duplicate symbolic graph entries are preserved.
 """
 function MPO_symbolic(
-  os::OpIDSum{N,C,Ti},
+  os::OpIDSum{N,W,Ti},
   sites::Vector{<:Index};
   basis_op_cache_vec=nothing,
   check_for_errors::Bool=true,
   combine_qn_sectors::Bool=false,
   output_level::Int=0,
   kwargs...,
-)::SymbolicMPO where {N,C,Ti}
-  _check_MPO_symbolic_kwargs(kwargs)
-
-  C <: Integer || throw(
-    ArgumentError("MPO_symbolic requires an OpIDSum with integer coefficient labels."),
-  )
-
-  internalize_symbolic_ids!(os)
-  prepare_opID_sum!(os, to_OpCacheVec(sites, basis_op_cache_vec); symbolic_coefficients=true)
+)::SymbolicMPO where {N,W<:AbstractWeight,Ti}
+  prepare_opID_sum!(os, to_OpCacheVec(sites, basis_op_cache_vec))
   check_for_errors && check_os_for_errors(os)
-  max_user_label = _max_user_label(os)
 
   label = "Constructing symbolic MPOGraph from $(length(os)) terms"
-  @time_if output_level 0 label g = MPOGraph(os; symbolic_coefficients=true)
+  @time_if output_level 0 label g = MPOGraph(os)
 
   llinks = Vector{Index}(undef, length(sites) + 1)
   if hasqns(sites)
@@ -654,7 +609,7 @@ function MPO_symbolic(
   end
 
   offsets = Vector{Vector{Int}}(undef, length(sites))
-  block_sparse_matrices = Vector{Vector{BlockSparseMatrix{SymbolicLocalMatrix{Ti}}}}(
+  block_sparse_matrices = Vector{Vector{BlockSparseMatrix{SymbolicLocalMatrix{W, Ti}}}}(
     undef, length(sites)
   )
 
@@ -671,12 +626,7 @@ function MPO_symbolic(
     output_level,
   )
 
-  return SymbolicMPO(offsets, block_sparse_matrices, sites, llinks, os.op_cache_vec, max_user_label)
-end
-
-function MPO_symbolic(os::OpSum, sites::Vector{<:Index}; kwargs...)
-  _check_MPO_symbolic_kwargs(kwargs)
-  throw(ArgumentError("MPO_symbolic supports OpIDSum input only; OpSum input is not supported."))
+  return SymbolicMPO(offsets, block_sparse_matrices, sites, llinks, os.op_cache_vec)
 end
 
 @doc """
diff --git a/src/OpIDSum.jl b/src/OpIDSum.jl
index 57a289e..8a8667b 100644
--- a/src/OpIDSum.jl
+++ b/src/OpIDSum.jl
@@ -190,7 +190,7 @@ mutable struct OpIDSum{N,C,Ti}
   op_cache_vec::OpCacheVec
   abs_tol::Float64
   modify!::FunctionWrappers.FunctionWrapper{
-    C,
+    Float64,
     Tuple{
       SubArray{
         OpID{Ti},
@@ -227,7 +227,7 @@ function OpIDSum{N,C,Ti}(
   data = Vector{NTuple{N,OpID{Ti}}}(undef, max_terms)
   terms = reinterpret(reshape, OpID{Ti}, data)
   wrapped_modify! = FunctionWrappers.FunctionWrapper{
-    C,
+    Float64,
     Tuple{
       SubArray{
         OpID{Ti},
@@ -263,7 +263,7 @@ Terms with `abs(scalar) <= abs_tol` are ignored when added.
 function OpIDSum{N,C,Ti}(
   max_terms::Int, op_cache_vec::OpCacheVec; abs_tol::Real=0
 )::OpIDSum{N,C,Ti} where {N,C,Ti}
-  return OpIDSum{N,C,Ti}(max_terms, op_cache_vec, ops -> 1; abs_tol)
+  return OpIDSum{N,C,Ti}(max_terms, op_cache_vec, ops -> 1.0; abs_tol)
 end
 
 """
@@ -306,7 +306,7 @@ term-modification callback associated with `os` is applied. Terms below
 merge duplicate terms or check capacity; callers must construct `os` with enough
 room for all accepted terms.
 """
-function ITensorMPS.add!(os::OpIDSum, scalar::Number, ops)::Nothing
+function ITensorMPS.add!(os::OpIDSum, scalar, ops)::Nothing
   abs(scalar) <= os.abs_tol && return nothing
 
   num_terms = Threads.atomic_add!(os.num_terms, 1) + 1
@@ -338,7 +338,7 @@ end
 
 Vararg convenience method for `add!`.
 """
-function ITensorMPS.add!(os::OpIDSum, scalar::Number, ops::OpID...)::Nothing
+function ITensorMPS.add!(os::OpIDSum, scalar, ops::OpID...)::Nothing
   return add!(os, scalar, ops)
 end
 
@@ -406,16 +406,8 @@ internal signed symbolic ids. In that mode, rewrite factors must be exactly
 converting the symbolic id into an ordinary numeric coefficient.
 """
 @timeit function rewrite_in_operator_basis!(
-  os::OpIDSum{N,C,Ti}, basis_op_cache_vec::OpCacheVec; symbolic_coefficients::Bool=false
+  os::OpIDSum{N,C,Ti}, basis_op_cache_vec::OpCacheVec
 )::Nothing where {N,C,Ti}
-  if symbolic_coefficients && !(C <: Integer)
-    throw(
-      ArgumentError(
-        "Symbolic operator-basis rewrites require integer internal coefficient ids.",
-      ),
-    )
-  end
-
   op_cache_vec = os.op_cache_vec
 
   function scale_by_first_nz!(matrix::Matrix)
@@ -471,11 +463,6 @@ converting the symbolic id into an ordinary numeric coefficient.
     end
   end
 
-  function apply_rewrite_factor(scalar, coeff)
-    symbolic_coefficients && return C(_scale_symbolic_weight(scalar, coeff))
-    return scalar * coeff
-  end
-
   Threads.@threads for i in eachindex(os)
     scalar, ops = os[i]
     term_is_zero = false
@@ -504,7 +491,7 @@ converting the symbolic id into an ordinary numeric coefficient.
         )
       end
 
-      scalar = apply_rewrite_factor(scalar, coeff)
+      scalar *= coeff
 
       ops[a] = OpID{Ti}(basis_id, ops[a].n)
       for k in (a + 1):b
@@ -654,11 +641,10 @@ scalars.
 """
 function prepare_opID_sum!(
   os::OpIDSum,
-  basis_op_cache_vec::Union{Nothing,OpCacheVec};
-  symbolic_coefficients::Bool=false,
+  basis_op_cache_vec::Union{Nothing,OpCacheVec}
 )::Nothing
   if !isnothing(basis_op_cache_vec)
-    rewrite_in_operator_basis!(os, basis_op_cache_vec; symbolic_coefficients)
+    rewrite_in_operator_basis!(os, basis_op_cache_vec)
   end
 
   return nothing
diff --git a/src/SymbolicLocalMatrix.jl b/src/SymbolicLocalMatrix.jl
index 150cca4..d1d5b15 100644
--- a/src/SymbolicLocalMatrix.jl
+++ b/src/SymbolicLocalMatrix.jl
@@ -10,153 +10,15 @@ the stored sign. The absolute value of `signed_local_op_id` is the local
 operator id; negative local operator ids indicate that a Jordan-Wigner string
 should be applied while evaluating the cached operator matrix.
 """
-const SymbolicLocalMatrix{Ti<:Integer} = Vector{Tuple{Int,Ti}}
-
-function _check_symbolic_weight_id(signed_weight_id::Integer)::Nothing
-  signed_weight_id == 0 && throw(
-    ArgumentError(
-      "Symbolic weight id 0 is invalid. Internal id 1 is reserved for the constant one.",
-    ),
-  )
-  return nothing
-end
-
-function _check_symbolic_local_op_id(signed_local_op_id::Integer)::Nothing
-  signed_local_op_id == 0 && throw(
-    ArgumentError("Symbolic local operator id 0 is invalid. Operator ids are one-based."),
-  )
-  return nothing
-end
-
-function _internal_symbolic_id(user_label::Integer)::Int
-  user_label > 0 ||
-    throw(ArgumentError("Symbolic coefficient labels must be greater than zero."))
-  return user_label + 1
-end
-
-"""
-    internalize_symbolic_ids!(os::OpIDSum{N,C,Ti}) where {N,C<:Integer,Ti}
-
-Map the integer scalar labels stored in `os` into the internal symbolic id
-space used by symbolic MPO construction.
-
-User label `k > 0` is stored as `k + 1`, reserving internal id `1` for the
-constant one. Labels less than or equal to zero are invalid. The conversion
-mutates `os` in place and does not alter its operator cache, tolerance, or
-modification callback.
-"""
-function internalize_symbolic_ids!(
-  os::OpIDSum{N,C,Ti}
-)::OpIDSum{N,C,Ti} where {N,C<:Integer,Ti}
-  for i in eachindex(os)
-    os.scalars[i] = C(_internal_symbolic_id(os.scalars[i]))
-  end
-  return os
-end
-
-function _max_user_label(signed_weight_id::Integer)::Int
-  _check_symbolic_weight_id(signed_weight_id)
-  return abs(signed_weight_id) - 1
-end
-
-function _max_user_label(terms::SymbolicLocalMatrix)::Int
-  max_user_label = 0
-  for (signed_weight_id, _) in terms
-    max_user_label = max(max_user_label, _max_user_label(signed_weight_id))
-  end
-  return max_user_label
-end
-
-function _max_user_label(os::OpIDSum)::Int
-  max_user_label = 0
-  for i in eachindex(os)
-    max_user_label = max(max_user_label, _max_user_label(os.scalars[i]))
-  end
-  return max_user_label
-end
+const SymbolicLocalMatrix{Weight<:AbstractWeight, Ti<:Integer} = Vector{Tuple{Weight, Ti}}
 
 function _append_symbolic_local_matrix_term!(
-  terms::SymbolicLocalMatrix{Ti},
-  signed_weight_id::Integer,
-  signed_local_op_id::Integer,
-)::SymbolicLocalMatrix{Ti} where {Ti}
-  _check_symbolic_weight_id(signed_weight_id)
-  _check_symbolic_local_op_id(signed_local_op_id)
-  push!(terms, (Int(signed_weight_id), Ti(signed_local_op_id)))
-  return terms
-end
-
-function _normalize_symbolic_local_matrix!(
-  terms::SymbolicLocalMatrix{Ti}
-)::SymbolicLocalMatrix{Ti} where {Ti}
-  for (signed_weight_id, signed_local_op_id) in terms
-    _check_symbolic_weight_id(signed_weight_id)
-    _check_symbolic_local_op_id(signed_local_op_id)
-  end
-
-  sort!(terms; by=term -> (term[2], abs(term[1]), term[1]))
-
-  read_index = 1
-  write_index = 1
-  while read_index <= length(terms)
-    signed_local_op_id = terms[read_index][2]
-    abs_weight_id = abs(terms[read_index][1])
-    num_positive = 0
-    num_negative = 0
-
-    while read_index <= length(terms) &&
-        terms[read_index][2] == signed_local_op_id &&
-        abs(terms[read_index][1]) == abs_weight_id
-      if terms[read_index][1] > 0
-        num_positive += 1
-      else
-        num_negative += 1
-      end
-      read_index += 1
-    end
-
-    multiplicity = num_positive - num_negative
-    signed_weight_id = multiplicity > 0 ? abs_weight_id : -abs_weight_id
-    for _ in 1:abs(multiplicity)
-      terms[write_index] = (signed_weight_id, signed_local_op_id)
-      write_index += 1
-    end
-  end
-
-  resize!(terms, write_index - 1)
-  return terms
-end
-
-function _scale_symbolic_weight(signed_weight_id::Integer, factor)::Int
-  _check_symbolic_weight_id(signed_weight_id)
-  isone(factor) && return Int(signed_weight_id)
-  factor == -one(factor) && return -Int(signed_weight_id)
-  throw(
-    ArgumentError(
-      "Unsupported symbolic rewrite factor $factor. " *
-      "Symbolic MPO construction supports only exact +1 and -1 factors.",
-    ),
-  )
-end
-
-function _weight_value(signed_weight_id::Integer, coefficients::AbstractVector)
-  abs_weight_id = abs(signed_weight_id)
-  abs_weight_id == 1 && return signed_weight_id > 0 ? 1.0 : -1.0
-
-  user_label = abs_weight_id - 1
-
-  value = coefficients[user_label]
-  return signed_weight_id > 0 ? value : -value
-end
-
-function _symbolic_local_matrix_eltype(
-  coefficients::AbstractVector, op_cache::Vector{OpInfo}
-)::Type
-  val_type = promote_type(Float64, eltype(coefficients))
-  for op_info in op_cache
-    val_type = promote_type(val_type, eltype(op_info.matrix))
-  end
-  return val_type
+  terms::SymbolicLocalMatrix{W, Ti},
+  weight::W,
+  signed_op_id::Ti,
+)::Nothing where {W, Ti}
+  push!(terms, (weight, signed_op_id))
+  return nothing
 end
 
 function _evaluate_symbolic_local_matrix!(
@@ -168,7 +30,7 @@ function _evaluate_symbolic_local_matrix!(
   for (signed_weight_id, signed_local_op_id) in terms
     local_op_id = abs(Int(signed_local_op_id))
     local_op = op_cache[local_op_id].matrix
-    weight = _weight_value(signed_weight_id, coefficients)
+    weight = substitute_weight(signed_weight_id, coefficients)
 
     needs_jw += signed_local_op_id < 0
     add_to_local_matrix!(
diff --git a/src/SymbolicWeight.jl b/src/SymbolicWeight.jl
new file mode 100644
index 0000000..ae1ac62
--- /dev/null
+++ b/src/SymbolicWeight.jl
@@ -0,0 +1,28 @@
+# TODO Rename file to AbstractWeight
+abstract type AbstractWeight end
+
+# Abs is used to drop terms in OpIDSum which have small coefficients.
+Base.abs(::AbstractWeight) = Inf64
+
+struct SimpleWeight <: AbstractWeight
+  id::Int
+end
+
+Base.one(::Type{SimpleWeight}) = SimpleWeight(typemax(Int))
+
+function Base.:*(s::SimpleWeight, x)::SimpleWeight
+  x == +1 && return s
+  x == -1 && return SimpleWeight(-s.id)
+  error("SimpleWeight can only be multiplied by ±1, not $x")
+end
+
+function substitute_weight(weight::SimpleWeight, coefficients::Vector{ValType})::ValType where {ValType}
+  abs_weight_id = abs(weight.id)
+
+  value = one(ValType)
+  if abs_weight_id != typemax(Int)
+    value = coefficients[abs_weight_id]
+  end
+
+  return weight.id > 0 ? value : -value
+end
\ No newline at end of file
diff --git a/src/large-graph-mpo-common.jl b/src/large-graph-mpo-common.jl
index 664baf9..591a60a 100644
--- a/src/large-graph-mpo-common.jl
+++ b/src/large-graph-mpo-common.jl
@@ -7,12 +7,16 @@ tensor storage.
 The outer vector is indexed by component-local `right_link`. Each inner
 `Dictionary` maps `left_link` to one local block value. Numeric construction
 uses dense local operator matrices as `BlockSparseMatrix{Matrix{C}}`, while the
-symbolic path stores term lists as `BlockSparseMatrix{SymbolicLocalMatrix{Ti}}`.
+symbolic path stores term lists as `BlockSparseMatrix{SymbolicLocalMatrix}`.
 `at_site!` later returns offsets that place component-local right-link ids into
 the full outgoing MPO bond.
 """
 BlockSparseMatrix{MatrixType} = Vector{Dictionary{Int,MatrixType}}
 
+function BlockSparseMatrix(::Type{MatrixType}, ncols::Int)::BlockSparseMatrix{MatrixType} where {MatrixType}
+  return [Dictionary{Int,MatrixType}() for _ in 1:ncols]
+end
+
 """
     MPOGraph{N,C,Ti}
 
@@ -219,13 +223,9 @@ identity incoming link and the operator emitted at site 1, while right vertices
 carry the remaining operator tuples.
 """
 @timeit function MPOGraph(
-  os::OpIDSum{N,C,Ti}; symbolic_coefficients::Bool=false
+  os::OpIDSum{N,C,Ti}
 )::MPOGraph{N,C,Ti} where {N,C,Ti}
   @assert size(os.terms, 1) == N
-  if symbolic_coefficients && !(C <: Integer)
-    throw(ArgumentError("Symbolic MPO graph construction requires integer coefficient ids."))
-  end
-
   ## Reverse the terms in the sum, ignoring trailing identity operators.
   Threads.@threads for i in 1:length(os)
     for j in N:-1:1
@@ -247,9 +247,8 @@ carry the remaining operator tuples.
   ## Combine duplicate numeric terms and remove terms which are below the tolerance.
   ## Symbolic ids are labels, not numeric weights, so preserve duplicate entries.
   nnz = 0
-  if symbolic_coefficients
+  if C <: AbstractWeight
     for i in eachindex(os)
-      _check_symbolic_weight_id(os.scalars[i])
       nnz += 1
       os.scalars[nnz] = os.scalars[i]
       os._data[nnz] = os._data[i]
@@ -410,7 +409,7 @@ known.
   alg::String;
   combine_qn_sectors::Bool,
   output_level::Int=0,
-)::Tuple{MPOGraph{N,C,Ti},Vector{Int},Vector{MatrixType},Index} where {MatrixType,N,C,Ti}
+)::Tuple{MPOGraph{N,C,Ti},Vector{Int},Vector{BlockSparseMatrix{MatrixType}},Index} where {MatrixType,N,C,Ti}
   has_qns = hasqns(sites)
 
   workspace = combine_duplicate_adjacent_right_vertices!(g, terms_eq_from(n + 1))
@@ -423,7 +422,7 @@ known.
   rank_of_cc = zeros(Int, nccs)
 
   ## The MPO tensor for each component.
-  matrix_of_cc = [_vertex_cover_matrix(MatrixType, 0) for _ in 1:nccs]
+  matrix_of_cc = [BlockSparseMatrix{MatrixType}() for _ in 1:nccs]
 
   ## The QN of each component
   qi_of_cc = Vector{Pair{QN,Int}}(undef, nccs)
diff --git a/src/large-graph-mpo-vertex-cover.jl b/src/large-graph-mpo-vertex-cover.jl
index fb56287..ffc7654 100644
--- a/src/large-graph-mpo-vertex-cover.jl
+++ b/src/large-graph-mpo-vertex-cover.jl
@@ -1,14 +1,3 @@
-function _vertex_cover_matrix(
-  ::Type{BlockSparseMatrix{Matrix{ValType}}}, rank::Int
-)::BlockSparseMatrix{Matrix{ValType}} where {ValType}
-  return [Dictionary{Int,Matrix{ValType}}() for _ in 1:rank]
-end
-
-function _vertex_cover_matrix(
-  ::Type{BlockSparseMatrix{SymbolicLocalMatrix{Ti}}}, rank::Int
-)::BlockSparseMatrix{SymbolicLocalMatrix{Ti}} where {Ti<:Integer}
-  return [Dictionary{Int,SymbolicLocalMatrix{Ti}}() for _ in 1:rank]
-end
 
 function _signed_symbolic_local_op_id(lv::LeftVertex, ::Type{Ti})::Ti where {Ti<:Integer}
   signed_local_op_id = lv.needs_JW_string ? -Int(lv.op_id) : Int(lv.op_id)
@@ -32,61 +21,22 @@ function _add_vertex_cover_term!(
 end
 
 function _add_vertex_cover_term!(
-  matrix::BlockSparseMatrix{SymbolicLocalMatrix{Ti}},
+  matrix::BlockSparseMatrix{SymbolicLocalMatrix{W, Ti}},
   m::Int,
   lv::LeftVertex,
-  signed_weight_id::Integer,
+  weight::W,
   op_cache::Vector{OpInfo},
   site_dim::Int,
-)::Nothing where {Ti<:Integer}
+)::Nothing where {W, Ti}
   terms = get!(matrix[m], lv.link) do
-    SymbolicLocalMatrix{Ti}()
+    SymbolicLocalMatrix{W, Ti}()
   end
   _append_symbolic_local_matrix_term!(
-    terms, signed_weight_id, _signed_symbolic_local_op_id(lv, Ti)
+    terms, weight, _signed_symbolic_local_op_id(lv, Ti)
   )
   return nothing
 end
 
-function _add_vertex_cover_unit_term!(
-  matrix::BlockSparseMatrix{Matrix{ValType}},
-  m::Int,
-  lv::LeftVertex,
-  op_cache::Vector{OpInfo},
-  site_dim::Int,
-  ::Type{C},
-)::Nothing where {ValType,C}
-  _add_vertex_cover_term!(matrix, m, lv, one(ValType), op_cache, site_dim)
-  return nothing
-end
-
-function _add_vertex_cover_unit_term!(
-  matrix::BlockSparseMatrix{SymbolicLocalMatrix{Ti}},
-  m::Int,
-  lv::LeftVertex,
-  op_cache::Vector{OpInfo},
-  site_dim::Int,
-  ::Type{C},
-)::Nothing where {Ti<:Integer,C<:Integer}
-  _add_vertex_cover_term!(matrix, m, lv, one(C), op_cache, site_dim)
-  return nothing
-end
-
-function _finalize_vertex_cover_matrix!(matrix::BlockSparseMatrix)::Nothing
-  return nothing
-end
-
-function _finalize_vertex_cover_matrix!(
-  matrix::BlockSparseMatrix{SymbolicLocalMatrix{Ti}}
-)::Nothing where {Ti}
-  for block in matrix
-    for terms in values(block)
-      _normalize_symbolic_local_matrix!(terms)
-    end
-  end
-  return nothing
-end
-
 """
     process_vertex_cover!(
       matrix_of_cc, rank_of_cc, next_edges_of_cc, g, ccs, n, sites, op_cache_vec
@@ -101,7 +51,7 @@ stored edge weights. On the final site, only the local tensor blocks are needed;
 otherwise this also builds `next_edges_of_cc` for the graph at site `n + 1`.
 """
 @timeit function process_vertex_cover!(
-  matrix_of_cc::Vector{MatrixType},
+  matrix_of_cc::Vector{BlockSparseMatrix{MatrixType}},
   rank_of_cc::Vector{Int},
   next_edges_of_cc::Vector{Matrix{Tuple{Vector{Int},Vector{C}}}},
   g::MPOGraph{N,C,Ti},
@@ -132,7 +82,7 @@ otherwise this also builds `next_edges_of_cc` for the graph at site `n + 1`.
     rank = length(left_cover) + length(right_cover)
     rank_of_cc[cc] = rank
 
-    matrix = _vertex_cover_matrix(MatrixType, rank)
+    matrix = BlockSparseMatrix(MatrixType, rank)
     matrix_of_cc[cc] = matrix
 
     ## Construct the tensor from the left cover.
@@ -140,7 +90,7 @@ otherwise this also builds `next_edges_of_cc` for the graph at site `n + 1`.
       lv_id = lvs_of_component[left_cover[m]]
       lv = left_vertex(g, lv_id)
 
-      _add_vertex_cover_unit_term!(matrix, m, lv, op_cache, site_dim, C)
+      _add_vertex_cover_term!(matrix, m, lv, one(C), op_cache, site_dim)
     end
 
     ## Construct the tensor from the right cover.
@@ -175,8 +125,6 @@ otherwise this also builds `next_edges_of_cc` for the graph at site `n + 1`.
       end
     end
 
-    _finalize_vertex_cover_matrix!(matrix)
-
     n == length(sites) && continue
 
     ## Preallocate space for next_edges
diff --git a/test/symbolic-local-matrix-tests.jl b/test/symbolic-local-matrix-tests.jl
index 72064ca..e1b6cd2 100644
--- a/test/symbolic-local-matrix-tests.jl
+++ b/test/symbolic-local-matrix-tests.jl
@@ -1,19 +1,15 @@
 using ITensorMPOConstruction
 using ITensorMPOConstruction:
   MPOGraph,
+  SimpleWeight,
   SymbolicLocalMatrix,
   _append_symbolic_local_matrix_term!,
   _evaluate_symbolic_local_matrix,
-  _internal_symbolic_id,
-  _max_user_label,
-  _normalize_symbolic_local_matrix!,
-  _scale_symbolic_weight,
-  _weight_value,
   combine_duplicate_adjacent_right_vertices!,
-  internalize_symbolic_ids!,
   left_size,
   left_vertex,
-  prepare_opID_sum!
+  prepare_opID_sum!,
+  substitute_weight
 using ITensorMPOConstruction: right_size, terms_eq_from
 using ITensors
 using ITensorMPS
@@ -22,42 +18,26 @@ using Test
 function test_symbolic_local_matrix_operations()::Nothing
   coefficients = [10.0, 20.0]
 
-  @test _weight_value(1, coefficients) == 1.0
-  @test _weight_value(-1, coefficients) == -1.0
-  @test _weight_value(3, coefficients) == coefficients[2]
-  @test _weight_value(-3, coefficients) == -coefficients[2]
+  @test substitute_weight(one(SimpleWeight), coefficients) == 1.0
+  @test substitute_weight(SimpleWeight(1), coefficients) == coefficients[1]
+  @test substitute_weight(SimpleWeight(-1), coefficients) == -coefficients[1]
+  @test substitute_weight(SimpleWeight(2), coefficients) == coefficients[2]
 
-  @test _internal_symbolic_id(1) == 2
-  @test _internal_symbolic_id(3) == 4
-  @test_throws ArgumentError _internal_symbolic_id(0)
-  @test_throws ArgumentError _internal_symbolic_id(-1)
+  @test SimpleWeight(3) * 1 == SimpleWeight(3)
+  @test SimpleWeight(3) * -1 == SimpleWeight(-3)
+  @test SimpleWeight(-3) * -1 == SimpleWeight(3)
+  @test_throws ErrorException SimpleWeight(3) * 2
+  @test_throws ErrorException SimpleWeight(3) * im
 
-  @test _scale_symbolic_weight(3, 1) == 3
-  @test _scale_symbolic_weight(3, -1) == -3
-  @test _scale_symbolic_weight(-3, -1) == 3
-  @test_throws ArgumentError _scale_symbolic_weight(3, 2)
-  @test_throws ArgumentError _scale_symbolic_weight(3, im)
-
-  terms = SymbolicLocalMatrix{Int}()
-  _append_symbolic_local_matrix_term!(terms, 3, -2)
-  @test terms == [(3, -2)]
-  @test _max_user_label(SymbolicLocalMatrix{Int}([(1, 1), (-3, 2), (5, -2)])) == 4
-
-  terms = SymbolicLocalMatrix{Int}([(3, 2), (-3, 2)])
-  _normalize_symbolic_local_matrix!(terms)
-  @test isempty(terms)
-
-  terms = SymbolicLocalMatrix{Int}([(3, 2), (3, 2)])
-  _normalize_symbolic_local_matrix!(terms)
-  @test terms == [(3, 2), (3, 2)]
-
-  terms = SymbolicLocalMatrix{Int}([(5, 3), (4, 2), (3, 2)])
-  _normalize_symbolic_local_matrix!(terms)
-  @test terms == [(3, 2), (4, 2), (5, 3)]
+  terms = SymbolicLocalMatrix{SimpleWeight,Int}()
+  _append_symbolic_local_matrix_term!(terms, SimpleWeight(1), -2)
+  @test terms == [(SimpleWeight(1), -2)]
 
   qubit_sites = siteinds("Qubit", 1)
   qubit_op_cache_vec = to_OpCacheVec(qubit_sites, [["I", "X"]])
-  duplicate_terms = SymbolicLocalMatrix{Int}([(3, 2), (3, 2)])
+  duplicate_terms = SymbolicLocalMatrix{SimpleWeight,Int}([
+    (SimpleWeight(2), 2), (SimpleWeight(2), 2)
+  ])
   @test _evaluate_symbolic_local_matrix(
     Float64, duplicate_terms, coefficients, qubit_op_cache_vec[1]
   ) ≈ 2 * coefficients[2] * qubit_op_cache_vec[1][2].matrix
@@ -65,10 +45,16 @@ function test_symbolic_local_matrix_operations()::Nothing
   fermion_sites = siteinds("Fermion", 1)
   fermion_op_cache_vec = to_OpCacheVec(fermion_sites, [["I", "C", "Cdag", "N"]])
   plain_matrix = _evaluate_symbolic_local_matrix(
-    Float64, SymbolicLocalMatrix{Int}([(3, 2)]), coefficients, fermion_op_cache_vec[1]
+    Float64,
+    SymbolicLocalMatrix{SimpleWeight,Int}([(SimpleWeight(2), 2)]),
+    coefficients,
+    fermion_op_cache_vec[1],
   )
   jw_matrix = _evaluate_symbolic_local_matrix(
-    Float64, SymbolicLocalMatrix{Int}([(3, -2)]), coefficients, fermion_op_cache_vec[1]
+    Float64,
+    SymbolicLocalMatrix{SimpleWeight,Int}([(SimpleWeight(2), -2)]),
+    coefficients,
+    fermion_op_cache_vec[1],
   )
 
   expected_plain_matrix = coefficients[2] * fermion_op_cache_vec[1][2].matrix
@@ -81,56 +67,6 @@ function test_symbolic_local_matrix_operations()::Nothing
   return nothing
 end
 
-function test_internalize_symbolic_ids()::Nothing
-  sites = siteinds("Qubit", 2)
-  op_cache_vec = to_OpCacheVec(sites, [["I", "X"], ["I", "Z"]])
-
-  X(n) = OpID(2, n)
-  Z(n) = OpID(2, n)
-
-  modify!(ops) = 1
-  os = OpIDSum{2,Int,Int}(3, op_cache_vec, modify!; abs_tol=0.25)
-  add!(os, 1, X(1))
-  add!(os, 2, Z(2))
-  add!(os, 3, X(1), Z(2))
-
-  original_op_cache_vec = os.op_cache_vec
-  original_abs_tol = os.abs_tol
-  original_modify! = os.modify!
-  original_scalars = os.scalars
-
-  @test internalize_symbolic_ids!(os) === os
-  @test original_scalars === os.scalars
-  @test os.scalars[1:3] == [2, 3, 4]
-  @test os.op_cache_vec === original_op_cache_vec
-  @test os.abs_tol == original_abs_tol
-  @test os.modify! === original_modify!
-
-  zero_label_os = OpIDSum{2,Int,Int}(1, op_cache_vec)
-  add!(zero_label_os, 1, X(1))
-  zero_label_os.scalars[1] = 0
-  @test_throws ArgumentError internalize_symbolic_ids!(zero_label_os)
-
-  negative_label_os = OpIDSum{2,Int,Int}(1, op_cache_vec)
-  add!(negative_label_os, -1, X(1))
-  @test_throws ArgumentError internalize_symbolic_ids!(negative_label_os)
-
-  return nothing
-end
-
-function test_negative_internal_symbolic_ids_remain_supported()::Nothing
-  coefficients = [7.0]
-
-  @test _weight_value(-2, coefficients) == -coefficients[1]
-
-  terms = SymbolicLocalMatrix{Int}()
-  _append_symbolic_local_matrix_term!(terms, -2, 2)
-  @test terms == [(-2, 2)]
-  @test _max_user_label(terms) == 1
-
-  return nothing
-end
-
 function symbolic_basis_rewrite_cache_vecs(factor)
   identity_matrix = [1.0 0.0; 0.0 1.0]
   x_matrix = [0.0 1.0; 1.0 0.0]
@@ -153,33 +89,31 @@ function symbolic_rewrite_opID_sum(factor, user_label::Int=7)
   scaled_identity(n) = OpID(2, n)
   X(n) = OpID(3, n)
 
-  os = OpIDSum{2,Int,Int}(1, op_cache_vec)
-  add!(os, user_label, scaled_identity(1), X(1))
-  internalize_symbolic_ids!(os)
+  os = OpIDSum{2,SimpleWeight,Int}(1, op_cache_vec)
+  add!(os, SimpleWeight(user_label), scaled_identity(1), X(1))
 
   return os, basis_op_cache_vec
 end
 
 function test_symbolic_basis_rewrite_sign_factors()::Nothing
   user_label = 7
-  internal_id = user_label + 1
 
   os, basis_op_cache_vec = symbolic_rewrite_opID_sum(-1.0, user_label)
-  prepare_opID_sum!(os, basis_op_cache_vec; symbolic_coefficients=true)
+  prepare_opID_sum!(os, basis_op_cache_vec)
 
   scalar, ops = os[1]
   nonzero_ops = [op for op in ops if op != zero(op)]
-  @test scalar == -internal_id
+  @test scalar == SimpleWeight(-user_label)
   @test nonzero_ops == [OpID(2, 1)]
   @test count(op -> op == zero(op), ops) == 1
   @test os.op_cache_vec === basis_op_cache_vec
 
   os, basis_op_cache_vec = symbolic_rewrite_opID_sum(1.0, user_label)
-  prepare_opID_sum!(os, basis_op_cache_vec; symbolic_coefficients=true)
+  prepare_opID_sum!(os, basis_op_cache_vec)
 
   scalar, ops = os[1]
   nonzero_ops = [op for op in ops if op != zero(op)]
-  @test scalar == internal_id
+  @test scalar == SimpleWeight(user_label)
   @test nonzero_ops == [OpID(2, 1)]
   @test count(op -> op == zero(op), ops) == 1
   @test os.op_cache_vec === basis_op_cache_vec
@@ -192,15 +126,14 @@ function test_symbolic_basis_rewrite_rejects_unsupported_factors()::Nothing
     os, basis_op_cache_vec = symbolic_rewrite_opID_sum(factor)
     err = nothing
     try
-      prepare_opID_sum!(os, basis_op_cache_vec; symbolic_coefficients=true)
+      prepare_opID_sum!(os, basis_op_cache_vec)
     catch e
       err = e
     end
 
     @test !isnothing(err)
     error_message = sprint(showerror, err)
-    @test occursin("ArgumentError", error_message)
-    @test occursin("Unsupported symbolic rewrite factor", error_message)
+    @test occursin("SimpleWeight can only be multiplied", error_message)
   end
 
   return nothing
@@ -212,15 +145,15 @@ function two_site_qubit_symbolic_fixture()
   return sites, op_cache_vec
 end
 
-function symbolic_terms_from_first_left_vertex(g)::SymbolicLocalMatrix{Int}
+function symbolic_terms_from_first_left_vertex(g)::SymbolicLocalMatrix{SimpleWeight,Int}
   @test left_size(g) == 1
   lv = left_vertex(g, 1)
   signed_local_op_id = lv.needs_JW_string ? -lv.op_id : lv.op_id
-  terms = SymbolicLocalMatrix{Int}()
+  terms = SymbolicLocalMatrix{SimpleWeight,Int}()
   for signed_weight_id in g.edge_weights_from_left[1]
-    _append_symbolic_local_matrix_term!(terms, signed_weight_id, signed_local_op_id)
+    _append_symbolic_local_matrix_term!(terms, signed_weight_id, Int(signed_local_op_id))
   end
-  return _normalize_symbolic_local_matrix!(terms)
+  return terms
 end
 
 function compact_duplicate_symbolic_right_vertices!(g)::Nothing
@@ -259,18 +192,17 @@ function test_symbolic_mpo_graph_preserves_duplicate_signed_ids()::Nothing
 
   X(n) = OpID(2, n)
 
-  os = OpIDSum{2,Int,Int}(2, op_cache_vec)
-  add!(os, 1, X(1), X(2))
-  add!(os, 2, X(1), X(2))
-  internalize_symbolic_ids!(os)
+  os = OpIDSum{2,SimpleWeight,Int}(2, op_cache_vec)
+  add!(os, SimpleWeight(1), X(1), X(2))
+  add!(os, SimpleWeight(2), X(1), X(2))
 
-  g = MPOGraph(os; symbolic_coefficients=true)
+  g = MPOGraph(os)
   @test length(g.edge_weights_from_left[1]) == 2
-  @test sort(g.edge_weights_from_left[1]) == [2, 3]
+  @test sort(getfield.(g.edge_weights_from_left[1], :id)) == [1, 2]
 
   compact_duplicate_symbolic_right_vertices!(g)
   terms = symbolic_terms_from_first_left_vertex(g)
-  @test terms == [(2, 2), (3, 2)]
+  @test terms == [(SimpleWeight(1), 2), (SimpleWeight(2), 2)]
 
   coefficients = [11.0, 17.0]
   evaluated = _evaluate_symbolic_local_matrix(Float64, terms, coefficients, op_cache_vec[1])
@@ -284,15 +216,14 @@ function test_symbolic_mpo_graph_preserves_duplicate_label_multiplicity()::Nothi
 
   X(n) = OpID(2, n)
 
-  os = OpIDSum{2,Int,Int}(2, op_cache_vec)
-  add!(os, 1, X(1), X(2))
-  add!(os, 1, X(1), X(2))
-  internalize_symbolic_ids!(os)
+  os = OpIDSum{2,SimpleWeight,Int}(2, op_cache_vec)
+  add!(os, SimpleWeight(1), X(1), X(2))
+  add!(os, SimpleWeight(1), X(1), X(2))
 
-  g = MPOGraph(os; symbolic_coefficients=true)
+  g = MPOGraph(os)
   compact_duplicate_symbolic_right_vertices!(g)
   terms = symbolic_terms_from_first_left_vertex(g)
-  @test terms == [(2, 2), (2, 2)]
+  @test terms == [(SimpleWeight(1), 2), (SimpleWeight(1), 2)]
 
   coefficient = 13.0
   evaluated = _evaluate_symbolic_local_matrix(Float64, terms, [coefficient], op_cache_vec[1])
@@ -327,16 +258,15 @@ function test_symbolic_mpo_graph_cancels_opposite_signed_ids()::Nothing
   X(n) = OpID(2, n)
   negX(n) = OpID(3, n)
 
-  os = OpIDSum{2,Int,Int}(2, op_cache_vec)
-  add!(os, 1, X(1), X(2))
-  add!(os, 1, negX(1), X(2))
-  internalize_symbolic_ids!(os)
-  prepare_opID_sum!(os, basis_op_cache_vec; symbolic_coefficients=true)
+  os = OpIDSum{2,SimpleWeight,Int}(2, op_cache_vec)
+  add!(os, SimpleWeight(1), X(1), X(2))
+  add!(os, SimpleWeight(1), negX(1), X(2))
+  prepare_opID_sum!(os, basis_op_cache_vec)
 
-  g = MPOGraph(os; symbolic_coefficients=true)
+  g = MPOGraph(os)
   compact_duplicate_symbolic_right_vertices!(g)
   terms = symbolic_terms_from_first_left_vertex(g)
-  @test isempty(terms)
+  @test terms == [(SimpleWeight(1), 2), (SimpleWeight(-1), 2)]
 
   evaluated = _evaluate_symbolic_local_matrix(Float64, terms, [7.0], basis_op_cache_vec[1])
   @test iszero(evaluated)
@@ -349,11 +279,6 @@ end
     test_symbolic_local_matrix_operations()
   end
 
-  @testset "internal symbolic ids" begin
-    test_internalize_symbolic_ids()
-    test_negative_internal_symbolic_ids_remain_supported()
-  end
-
   @testset "basis rewrite" begin
     test_symbolic_basis_rewrite_sign_factors()
     test_symbolic_basis_rewrite_rejects_unsupported_factors()
diff --git a/test/symbolic-mpo-tests.jl b/test/symbolic-mpo-tests.jl
index 6a106e8..534bd55 100644
--- a/test/symbolic-mpo-tests.jl
+++ b/test/symbolic-mpo-tests.jl
@@ -1,5 +1,5 @@
 using ITensorMPOConstruction
-using ITensorMPOConstruction: MPO_symbolic, SymbolicMPO
+using ITensorMPOConstruction: MPO_symbolic, SimpleWeight, SymbolicMPO
 using ITensors
 using ITensorMPS
 using Test
@@ -17,10 +17,10 @@ function transverse_field_ising_symbolic_fixture()
   X(n) = OpID(2, n)
   Z(n) = OpID(3, n)
 
-  os = OpIDSum{2,Int,Int}(3, op_cache_vec)
-  add!(os, 1, X(1), X(2))
-  add!(os, 2, Z(1))
-  add!(os, 3, Z(2))
+  os = OpIDSum{2,SimpleWeight,Int}(3, op_cache_vec)
+  add!(os, SimpleWeight(1), X(1), X(2))
+  add!(os, SimpleWeight(2), Z(1))
+  add!(os, SimpleWeight(3), Z(2))
 
   return os, sites
 end
@@ -45,10 +45,10 @@ function qn_number_symbolic_fixture()
 
   Nop(n) = OpID(2, n)
 
-  os = OpIDSum{2,Int,Int}(3, op_cache_vec)
-  add!(os, 1, Nop(1))
-  add!(os, 2, Nop(2))
-  add!(os, 3, Nop(1), Nop(2))
+  os = OpIDSum{2,SimpleWeight,Int}(3, op_cache_vec)
+  add!(os, SimpleWeight(1), Nop(1))
+  add!(os, SimpleWeight(2), Nop(2))
+  add!(os, SimpleWeight(3), Nop(1), Nop(2))
 
   return os, sites
 end
@@ -79,8 +79,8 @@ function test_mpos_approx_equal(A::MPO, B::MPO; tol::Real=1e-10)::Nothing
   return nothing
 end
 
-function symbolic_mpo_terms(sym::SymbolicMPO)::Vector{Tuple{Int,Int}}
-  terms = Tuple{Int,Int}[]
+function symbolic_mpo_terms(sym::SymbolicMPO)::Vector
+  terms = []
   for site_matrices in sym.block_sparse_matrices
     for matrix_of_component in site_matrices
       for block in matrix_of_component
@@ -102,28 +102,12 @@ function test_symbolic_mpo_construction_metadata()::Nothing
   @test length(sym.offsets) == length(sites)
   @test length(sym.block_sparse_matrices) == length(sites)
   @test length(sym.llinks) == length(sites) + 1
-  @test sym.max_user_label == 3
   @test sym.op_cache_vec === os.op_cache_vec
   @test !isempty(symbolic_mpo_terms(sym))
 
   return nothing
 end
 
-function test_symbolic_mpo_rejects_unsupported_public_inputs()::Nothing
-  os, sites = transverse_field_ising_symbolic_fixture()
-  @test_throws ArgumentError MPO_symbolic(os, sites; alg="QR")
-
-  noninteger_os = OpIDSum{2,Float64,Int}(1, os.op_cache_vec)
-  add!(noninteger_os, 1.0, OpID(2, 1), OpID(2, 2))
-  @test_throws ArgumentError MPO_symbolic(noninteger_os, sites)
-
-  op_sum = OpSum()
-  op_sum += "X", 1
-  @test_throws ArgumentError MPO_symbolic(op_sum, sites)
-
-  return nothing
-end
-
 function test_symbolic_mpo_fresh_instantiation_matches_numeric()::Nothing
   os, sites = transverse_field_ising_symbolic_fixture()
   sym = MPO_symbolic(os, sites)
@@ -149,9 +133,9 @@ function test_symbolic_mpo_duplicate_labels_instantiate_like_numeric()::Nothing
 
   X(n) = OpID(2, n)
 
-  symbolic_os = OpIDSum{2,Int,Int}(2, op_cache_vec)
-  add!(symbolic_os, 1, X(1), X(2))
-  add!(symbolic_os, 2, X(1), X(2))
+  symbolic_os = OpIDSum{2,SimpleWeight,Int}(2, op_cache_vec)
+  add!(symbolic_os, SimpleWeight(1), X(1), X(2))
+  add!(symbolic_os, SimpleWeight(2), X(1), X(2))
 
   sym = MPO_symbolic(symbolic_os, sites)
   coefficients = [2.0, 3.0]
@@ -189,7 +173,7 @@ function test_symbolic_mpo_qn_fresh_instantiation()::Nothing
   zero_mpo = instantiate_MPO(sym, [0.0, 0.0, 0.0]; splitblocks=true)
 
   @test linkdims(zero_mpo) == linkdims(nonzero_mpo)
-  @test_throws ArgumentError instantiate_MPO(sym, [1.0, 2.0])
+  @test_throws BoundsError instantiate_MPO(sym, [1.0, 2.0])
 
   return nothing
 end
@@ -200,11 +184,7 @@ function test_symbolic_mpo_missing_coefficients_rejected()::Nothing
   coefficients = [1.0, 2.0, 3.0]
   short_coefficients = [1.0, 2.0]
 
-  H = instantiate_MPO(sym, coefficients; splitblocks=true)
-
-  @test_throws ArgumentError instantiate_MPO(sym, short_coefficients; splitblocks=true)
-  @test_throws ArgumentError instantiate_MPO!(H, sym, short_coefficients)
-  @test_throws ArgumentError instantiate_MPO(H, sym, short_coefficients)
+  @test_throws BoundsError instantiate_MPO(sym, short_coefficients; splitblocks=true)
 
   return nothing
 end
@@ -284,8 +264,8 @@ function test_symbolic_mpo_uses_negative_local_op_id_for_jw_terms()::Nothing
   C(n) = OpID(2, n)
   Cdag(n) = OpID(3, n)
 
-  os = OpIDSum{2,Int,Int}(1, op_cache_vec)
-  add!(os, 1, Cdag(1), C(2))
+  os = OpIDSum{2,SimpleWeight,Int}(1, op_cache_vec)
+  add!(os, SimpleWeight(1), Cdag(1), C(2))
 
   sym = MPO_symbolic(os, sites)
   terms = symbolic_mpo_terms(sym)
@@ -304,7 +284,6 @@ end
   end
 
   @testset "public API rejections" begin
-    test_symbolic_mpo_rejects_unsupported_public_inputs()
     test_symbolic_mpo_missing_coefficients_rejected()
     test_symbolic_mpo_incompatible_template_rejected()
   end