diff --git a/Project.toml b/Project.toml
index 690c757..5656970 100644
--- a/Project.toml
+++ b/Project.toml
@@ -18,7 +18,7 @@ TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 [compat]
 Dictionaries = "0.4.6"
 FunctionWrappers = "1.1.3"
-ITensorMPS = "0.1, 0.2, 0.3"
+ITensorMPS = "0.1, 0.2, 0.3, 0.4"
 ITensors = "0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9"
 LinearAlgebra = "1.9, 1.10"
 Printf = "1"
diff --git a/docs/src/documentation/MPO_new.md b/docs/src/documentation/MPO_new.md
index 49ddafc..35cb60b 100644
--- a/docs/src/documentation/MPO_new.md
+++ b/docs/src/documentation/MPO_new.md
@@ -3,7 +3,7 @@ The main exported function is `MPO_new`, which takes an `OpSum` and transforms i
 
 1. The operator must be expressed as a sum of products of single-site operators. For example, a CNOT could not appear in the sum since it is a two-site operator.
 
-2. When dealing with fermionic systems, the parity of each term in the sum must be even. That is, the combined number of creation and annihilation operators in each term must be even. It should be possible to relax this constraint.
+2. When dealing with fermionic systems, every term in the sum must have the same fermion parity. The common parity may be even or odd, but even- and odd-parity terms cannot be mixed in a single sum.
 
 3. Each term in the sum-of-products representation can only have a single operator acting on any one site. For example, a term such as ``\mathbf{X}^{(1)} \mathbf{X}^{(1)}`` is not allowed. However, there is a preprocessing utility that can automatically replace ``\mathbf{X}^{(1)} \mathbf{X}^{(1)}`` with ``\mathbf{I}^{(1)}``.
 
diff --git a/examples/electronic-structure.jl b/examples/electronic-structure.jl
index 876af39..ec1d455 100644
--- a/examples/electronic-structure.jl
+++ b/examples/electronic-structure.jl
@@ -145,10 +145,9 @@ for alg in ("VC", "QR")
       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)
diff --git a/src/MPOConstruction.jl b/src/MPOConstruction.jl
index 5a8edd8..87c2445 100644
--- a/src/MPOConstruction.jl
+++ b/src/MPOConstruction.jl
@@ -166,19 +166,10 @@ function MPO_new(
   basis_op_cache_vec=nothing,
   check_for_errors::Bool=true,
   checkflux::Bool=true,
-  splitblocks::Union{Bool,Nothing}=nothing,
+  splitblocks::Bool=true,
   output_level::Int=0,
   kwargs...,
 )::MPO
-  if isnothing(splitblocks) # TODO: Remove warning some time after v0.2.1 release
-    splitblocks = true
-    hasqns(sites) && Base.depwarn(
-      "`splitblocks` not specified. The default is `true`, which is a change from prior behavior.",
-      :MPO_new;
-      force=true,
-    )
-  end
-
   prepare_opID_sum!(os, to_OpCacheVec(sites, basis_op_cache_vec))
   check_for_errors && check_os_for_errors(os)
 
diff --git a/src/OpIDSum.jl b/src/OpIDSum.jl
index e02eef5..a02d926 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
 
@@ -578,22 +576,25 @@ Validate structural assumptions for `os`.
 
 This checks that the cached local operators are linearly independent on every
 site and that, within every term, operators are sorted by site, at most one
-operator acts on each site, all terms carry the same total QN flux, and every
-term has even fermion parity.
+operator acts on each site, and all terms carry the same total QN flux and
+fermion parity. The common fermion parity may be even or odd.
 """
 @timeit function check_os_for_errors(os::OpIDSum)::Nothing
   flux_of_first_term = nothing
+  parity_of_first_term = nothing
   for i in eachindex(os)
     _, ops = os[i]
 
     flux = QN()
-    fermion_parity = 0
+    fermion_parity = false
     for j in eachindex(ops)
       opj = ops[j]
       opj == zero(opj) && continue
 
       flux += os.op_cache_vec[opj.n][opj.id].qn_flux
-      fermion_parity += os.op_cache_vec[opj.n][opj.id].is_fermionic
+      fermion_parity = xor(
+        fermion_parity, os.op_cache_vec[opj.n][opj.id].is_fermionic
+      )
 
       if j < length(ops)
         ops[j + 1] == zero(ops[j + 1]) && continue
@@ -613,7 +614,19 @@ term has even fermion parity.
       )
     end
 
-    mod(fermion_parity, 2) != 0 && error("Odd parity fermion terms not supported: $ops")
+    if isnothing(parity_of_first_term)
+      parity_of_first_term = fermion_parity
+    else
+      if fermion_parity != parity_of_first_term
+        first_parity_name = parity_of_first_term ? "odd" : "even"
+        current_parity_name = fermion_parity ? "odd" : "even"
+        error(
+          "Inconsistent fermion parity found!\n" *
+          "  Term 1 has $first_parity_name fermion parity\n" *
+          "  Term $i has $current_parity_name fermion parity",
+        )
+      end
+    end
   end
 
   check_op_cache_vec_linearly_independent(os.op_cache_vec)
diff --git a/test/Project.toml b/test/Project.toml
index 09c7811..a359852 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -10,7 +10,7 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 [compat]
 Dictionaries = "0.4.6"
 Graphs = "1.14.0"
-ITensorMPS = "0.1, 0.2, 0.3"
+ITensorMPS = "0.1, 0.2, 0.3, 0.4"
 ITensors = "0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9"
 LinearAlgebra = "1.9, 1.10"
 Random = "1"
diff --git a/test/mponew_h1_accuracy.jl b/test/mponew_h1_accuracy.jl
index 10c1a2a..3eeb15e 100644
--- a/test/mponew_h1_accuracy.jl
+++ b/test/mponew_h1_accuracy.jl
@@ -125,7 +125,7 @@ end
   @test maximum(abs, H .- transpose(H)) <= 1e-12
 
   sites = siteinds("Fermion", N; conserve_qns=true, conserve_nf=true)
-  H_new = MPO_new(build_h1_opsum(terms), sites; tol=1.0, splitblocks=true) # TODO: remove splitblocks=true after removing warning
+  H_new = MPO_new(build_h1_opsum(terms), sites; tol=1.0)
 
   tests = make_test_orbitals(N; n_basis=8, n_random=4, seed=1234)
   max_delta = 0.0
diff --git a/test/ops-tests.jl b/test/ops-tests.jl
index ac6f70d..501d634 100644
--- a/test/ops-tests.jl
+++ b/test/ops-tests.jl
@@ -187,6 +187,28 @@ function test_check_os_for_errors_rejects_linearly_dependent_cache()::Nothing
   return nothing
 end
 
+function test_check_os_for_errors_fermion_parity()::Nothing
+  sites = siteinds("Fermion", 3)
+  op_cache_vec =
+    to_OpCacheVec(sites, [["I", "C", "Cdag", "N"] for _ in eachindex(sites)])
+
+  C(n) = OpID(2, n)
+  Cdag(n) = OpID(3, n)
+  N(n) = OpID(4, n)
+
+  odd_os = OpIDSum{3,Float64,Int}(2, op_cache_vec)
+  add!(odd_os, 1.0, C(1))
+  add!(odd_os, 1.0, Cdag(1), C(2), C(3))
+  @test isnothing(check_os_for_errors(odd_os))
+
+  mixed_os = OpIDSum{2,Float64,Int}(2, op_cache_vec)
+  add!(mixed_os, 1.0, C(1))
+  add!(mixed_os, 1.0, N(2))
+  @test_throws "Inconsistent fermion parity found!" check_os_for_errors(mixed_os)
+
+  return nothing
+end
+
 @testset "Ops" begin
   test_are_equal()
 
@@ -201,4 +223,6 @@ end
   test_rewrite_in_operator_basis_zero()
 
   test_check_os_for_errors_rejects_linearly_dependent_cache()
+
+  test_check_os_for_errors_fermion_parity()
 end
diff --git a/test/test-MPOConstruction.jl b/test/test-MPOConstruction.jl
index 12a8212..43d4a00 100644
--- a/test/test-MPOConstruction.jl
+++ b/test/test-MPOConstruction.jl
@@ -254,7 +254,7 @@ function test_non_zero_flux(alg::String)::Nothing
   os = OpSum()
   os .+= 1.0, "S-", 1, "S-", 4
   os .+= 1.0, "S-", 2, "S-", 3
-  O = MPO_new(os, sites; alg, splitblocks=true) # TODO: remove splitblocks = true after warning.
+  O = MPO_new(os, sites; alg)
 
   @test flux(O) == QN(("Number", 2))
 
@@ -265,7 +265,7 @@ function test_non_zero_flux(alg::String)::Nothing
     os = OpSum()
     os .+= 1.0, "S-", 1
     os .+= 1.0, "S-", 2, "S-", 3
-    @test_throws "Inconsistent flux found!" MPO_new(os, sites; alg, splitblocks=true) # TODO: remove splitblocks = true after warning.
+    @test_throws "Inconsistent flux found!" MPO_new(os, sites; alg)
   end
 
   return nothing
@@ -348,6 +348,36 @@ function test_Fermi_Hubbard(
   return nothing
 end
 
+function test_odd_fermion_parity(alg::String, conserve_qns::Bool)::Nothing
+  sites = siteinds("Fermion", 5; conserve_qns)
+
+  os = OpSum{Float64}()
+  os .+= 0.7, "Cdag", 1
+  os .+= -1.2, "Cdag", 3
+  os .+= 0.4, "Cdag", 5
+  os .+= 0.9, "Cdag", 1, "Cdag", 3, "C", 5
+
+  mpo = MPO_new(os, sites; alg)
+
+  I(n) = op("I", sites[n])
+  F(n) = op("F", sites[n])
+  C(n) = op("C", sites[n])
+  Cdag(n) = op("Cdag", sites[n])
+
+  exact = 0.7 * Cdag(1) * I(2) * I(3) * I(4) * I(5)
+  exact += -1.2 * F(1) * F(2) * Cdag(3) * I(4) * I(5)
+  exact += 0.4 * F(1) * F(2) * F(3) * F(4) * Cdag(5)
+  exact += 0.9 * Cdag(1) * I(2) * op("Cdag * F", sites[3]) * F(4) * C(5)
+
+  @test norm(prod(mpo) - exact) < 1e-12
+  if conserve_qns
+    @test flux(mpo) == flux(Cdag(1))
+    @test !iszero(flux(mpo))
+  end
+
+  return nothing
+end
+
 ITensors.op(::OpName"00", ::SiteType"Qubit") = [1 0; 0 0]
 
 ITensors.op(::OpName"01", ::SiteType"Qubit") = [0 1; 0 0]
@@ -409,6 +439,11 @@ end
 
     @testset "$alg: non zero flux" test_non_zero_flux(alg)
 
+    @testset "$alg: odd fermion parity" begin
+      test_odd_fermion_parity(alg, false)
+      test_odd_fermion_parity(alg, true)
+    end
+
     @testset "$alg: qft" begin
       test_qft(6, false, alg)
       test_qft(6, true, alg)