diff --git a/Project.toml b/Project.toml
index 57d2753c..bb2537d3 100644
--- a/Project.toml
+++ b/Project.toml
@@ -24,6 +24,7 @@ NamedGraphs = "678767b0-92e7-4007-89e4-4527a8725b19"
 PackageExtensionCompat = "65ce6f38-6b18-4e1d-a461-8949797d7930"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
 SerializedElementArrays = "d3ce8812-9567-47e9-a7b5-65a6d70a3065"
 SimpleTraits = "699a6c99-e7fa-54fc-8d76-47d257e15c1d"
 SparseArrayKit = "a9a3c162-d163-4c15-8926-b8794fbefed2"
diff --git a/src/beliefpropagationdmrg/bp_dmrg.jl b/src/beliefpropagationdmrg/bp_dmrg.jl
new file mode 100644
index 00000000..bc02dc10
--- /dev/null
+++ b/src/beliefpropagationdmrg/bp_dmrg.jl
@@ -0,0 +1,55 @@
+using NamedGraphs.GraphsExtensions: is_tree
+using NamedGraphs.PartitionedGraphs: partitionvertices, partitionedges, PartitionEdge
+using ITensorNetworks: ITensorNetwork, QuadraticFormNetwork, BeliefPropagationCache, update, default_message_update, delete_messages
+using ITensors: scalar
+
+include("utils.jl")
+include("bp_extracter.jl")
+include("bp_inserter.jl")
+include("bp_updater.jl")
+include("graphsextensions.jl")
+
+default_bp_update_kwargs(ψ::ITensorNetwork) = is_tree(ψ) ? (;) : (; maxiter = 50, tol = 1e-8)
+
+function initialize_cache(ψ_init::ITensorNetwork; cache_update_kwargs = default_bp_update_kwargs(ψ_init))
+    ψ = copy(ψ_init)
+    ψIψ = QuadraticFormNetwork(ψ)
+    ψIψ_bpc = BeliefPropagationCache(ψIψ)
+    return (ψ, ψIψ_bpc)
+end
+
+function bp_dmrg(ψ_init::ITensorNetwork, H::OpSum; nsites = 1, no_sweeps = 1, bp_update_kwargs = default_bp_update_kwargs(ψ_init))
+
+    L = length(vertices(ψ_init))
+    state, ψIψ_bpc = initialize_cache(ψ_init)
+    state_vertices, state_edges = collect(vertices(state)), edges(state)
+    if nsites == 1
+        regions = [[v] for v in vcat(state_vertices, reverse(state_vertices))]
+    else 
+        regions = [[src(e), dst(e)] for e in vcat(state_edges, reverse(state_edges))]
+    end
+    state, ψIψ_bpc = renormalize_update_norm_cache(state, ψIψ_bpc; cache_update_kwargs = bp_update_kwargs)
+
+    energy = sum(expect(state, H; alg="bp", (cache!)=Ref(ψIψ_bpc))) / L
+    println("Initial energy density is $(energy)")
+    energies = [energy]
+
+    for i in 1:no_sweeps
+        println("Beginning sweep $i")
+        for region in regions
+            println("Updating vertex $region")
+
+            local_state, ∂ψOψ_bpc_∂rs, sqrt_mts, inv_sqrt_mts, ψ, ψIψ_bpc = bp_extracter_V3(state, H, ψIψ_bpc, region)
+
+            local_state, energy = bp_eigsolve_updater(local_state, ∂ψOψ_bpc_∂rs, sqrt_mts, inv_sqrt_mts)
+
+            state, ψIψ_bpc = bp_inserter(state, ψIψ_bpc, local_state, region; bp_update_kwargs)
+
+            energy = sum(expect(state, H; alg="bp", (cache!)=Ref(ψIψ_bpc))) / L
+            append!(energies, energy)
+            println("Current energy density is $(energy)")
+        end
+    end
+
+    return state, energies
+end
diff --git a/src/beliefpropagationdmrg/bp_extracter.jl b/src/beliefpropagationdmrg/bp_extracter.jl
new file mode 100644
index 00000000..00b6f099
--- /dev/null
+++ b/src/beliefpropagationdmrg/bp_extracter.jl
@@ -0,0 +1,83 @@
+using ITensors: scalartype, which_op, name, names
+using ITensorNetworks:
+  ket_vertices, bra_vertices, tensornetwork, default_message_update, operator_network
+using ITensorNetworks.ITensorsExtensions: map_eigvals
+
+function effective_environments(state::ITensorNetwork, H::OpSum, ψIψ_bpc::BeliefPropagationCache, region)
+  environments = Vector{ITensor}[]
+  s = siteinds(state)
+  op_tensors = Vector{ITensor}(H, s)
+  for (i, term) in enumerate(H)
+    if length(sites(term)) == 1 && !iszero(first(term.args))
+      path = a_star(state, only(sites(term)), only(region))
+      ψOψ_qf = QuadraticFormNetwork(state)
+
+      ψOψ_qf[(only(sites(term)), "operator")] = op_tensors[i]
+    elseif length(sites(term)) == 2 && !iszero(first(term.args))
+      v1, v2 = first(sites(term)), last(sites(term))
+      path = a_star(state, v1, only(region))
+      if v2 ∉ src.(path) && v2 ∉ dst.(path)
+        prepend!(path, NamedEdge[NamedEdge(v2 => v1)])
+      end
+      Ov1, Ov2 = factorize_svd(op_tensors[i], s[v1], s[v1]'; cutoff = 1e-16)
+      ψOψ_qf = QuadraticFormNetwork(state)
+      ψOψ_qf[(v1, "operator")] = Ov1
+      ψOψ_qf[(v2, "operator")] = Ov2
+    else 
+      path = nothing
+    end
+    if !isnothing(path)
+      env = ITensor(1.0)
+      path = vcat(src.(path))
+      for v in path
+        env *= get_local_term(ψOψ_qf, v)
+        vns = neighbors(state, v)
+        for vn in vns
+          if vn ∉ path && vn != only(region)
+            env *= only(message(ψIψ_bpc, PartitionEdge(vn => v)))
+          end
+        end
+      end
+      for vn in neighbors(state, only(region))
+        if vn ∉ path
+          env *= only(message(ψIψ_bpc, PartitionEdge(vn => only(region))))
+        end
+      end
+      env *= ψOψ_qf[(only(region), "operator")]
+
+      push!(environments, [env])
+    end
+  end
+
+  return environments
+end
+
+function bp_extracter(
+  ψ::AbstractITensorNetwork,
+  H::OpSum,
+  ψIψ_bpc::BeliefPropagationCache,
+  region;
+  regularization=10 * eps(scalartype(ψ)),
+  ishermitian=true,
+)
+
+  form_network = tensornetwork(ψIψ_bpc)
+  form_ket_vertices, form_bra_vertices = ket_vertices(form_network, region),
+  bra_vertices(form_network, region)
+
+  ∂ψOψ_bpc_∂rs = effective_environments(ψ, H, ψIψ_bpc, region)
+  state = prod([ψ[v] for v in region])
+  messages = environment(ψIψ_bpc, partitionvertices(ψIψ_bpc, form_ket_vertices))
+  f_sqrt = sqrt ∘ (x -> x + regularization)
+  f_inv_sqrt = inv ∘ sqrt ∘ (x -> x + regularization)
+  sqrt_mts =
+    map_eigvals.(
+      (f_sqrt,), messages, first.(inds.(messages)), last.(inds.(messages)); ishermitian
+    )
+  inv_sqrt_mts =
+    map_eigvals.(
+      (f_inv_sqrt,), messages, first.(inds.(messages)), last.(inds.(messages)); ishermitian
+    )
+
+  return state, ∂ψOψ_bpc_∂rs, sqrt_mts, inv_sqrt_mts, ψ, ψIψ_bpc
+end
\ No newline at end of file
diff --git a/src/beliefpropagationdmrg/bp_inserter.jl b/src/beliefpropagationdmrg/bp_inserter.jl
new file mode 100644
index 00000000..70aa91ad
--- /dev/null
+++ b/src/beliefpropagationdmrg/bp_inserter.jl
@@ -0,0 +1,101 @@
+using ITensorNetworks: update_factors, edge_tag
+using ITensors: uniqueinds, factorize_svd, factorize
+using LinearAlgebra: norm
+
+function renormalize_update_norm_cache(
+  ψ::ITensorNetwork,
+  ψOψ_bpcs::Vector{<:BeliefPropagationCache},
+  ψIψ_bpc::BeliefPropagationCache;
+  cache_update_kwargs,
+)
+  ψ = copy(ψ)
+  ψIψ_bpc = delete_messages(ψIψ_bpc)
+  ψOψ_bpcs = delete_messages.(ψOψ_bpcs)
+  ψIψ_bpc = update(ψIψ_bpc; cache_update_kwargs...)
+  ψIψ_bpc = renormalize_messages(ψIψ_bpc)
+  qf = tensornetwork(ψIψ_bpc)
+
+  for v in vertices(ψ)
+    v_ket, v_bra = ket_vertex(qf, v), bra_vertex(qf, v)
+    pv = only(partitionvertices(ψIψ_bpc, [v_ket]))
+    vn = region_scalar(ψIψ_bpc, pv)
+    state = (1.0 / sqrt(vn)) * ψ[v]
+    state_dag = copy(dag(state))
+    state_dag = replaceinds(
+      state_dag, inds(state_dag), dual_index_map(qf).(inds(state_dag))
+    )
+    vertices_states = Dictionary([v_ket, v_bra], [state, state_dag])
+    ψOψ_bpcs = update_factors.(ψOψ_bpcs, (vertices_states,))
+    ψIψ_bpc = update_factors(ψIψ_bpc, vertices_states)
+    ψ[v] = state
+  end
+
+  return ψ, ψOψ_bpcs, ψIψ_bpc
+end
+
+function renormalize_update_norm_cache(
+    ψ::ITensorNetwork,
+    ψIψ_bpc::BeliefPropagationCache;
+    cache_update_kwargs,
+  )
+    ψ = copy(ψ)
+    ψIψ_bpc = delete_messages(ψIψ_bpc)
+    ψIψ_bpc = update(ψIψ_bpc; cache_update_kwargs...)
+    ψIψ_bpc = renormalize_messages(ψIψ_bpc)
+    qf = tensornetwork(ψIψ_bpc)
+  
+    for v in vertices(ψ)
+      v_ket, v_bra = ket_vertex(qf, v), bra_vertex(qf, v)
+      pv = only(partitionvertices(ψIψ_bpc, [v_ket]))
+      vn = region_scalar(ψIψ_bpc, pv)
+      state = (1.0 / sqrt(vn)) * ψ[v]
+      state_dag = copy(dag(state))
+      state_dag = replaceinds(
+        state_dag, inds(state_dag), dual_index_map(qf).(inds(state_dag))
+      )
+      vertices_states = Dictionary([v_ket, v_bra], [state, state_dag])
+      ψIψ_bpc = update_factors(ψIψ_bpc, vertices_states)
+      ψ[v] = state
+    end
+  
+    return ψ, ψIψ_bpc
+end
+
+#TODO: Add support for nsites = 2
+function bp_inserter(
+    ψ::AbstractITensorNetwork,
+    ψIψ_bpc::BeliefPropagationCache,
+    state::ITensor,
+    region;
+    nsites::Int64=1,
+    bp_update_kwargs,
+    kwargs...,
+  )
+    ψ = copy(ψ)
+    form_network = tensornetwork(ψIψ_bpc)
+    if length(region) == 1
+      states = ITensor[state]
+    else
+      error("Region lengths of more than 1 not supported for now")
+    end
+  
+    for (state, v) in zip(states, region)
+      ψ[v] = state
+      state_dag = copy(ψ[v])
+      state_dag = replaceinds(
+        dag(state_dag), inds(state_dag), dual_index_map(form_network).(inds(state_dag))
+      )
+      form_bra_v, form_op_v, form_ket_v = bra_vertex(form_network, v),
+      operator_vertex(form_network, v),
+      ket_vertex(form_network, v)
+      vertices_states = Dictionary([form_ket_v, form_bra_v], [state, state_dag])
+      ψIψ_bpc = update_factors(ψIψ_bpc, vertices_states)
+    end
+  
+    ψ, ψIψ_bpc = renormalize_update_norm_cache(
+      ψ, ψIψ_bpc; cache_update_kwargs=bp_update_kwargs
+    )
+  
+    return ψ, ψIψ_bpc
+  end
+  
diff --git a/src/beliefpropagationdmrg/bp_updater.jl b/src/beliefpropagationdmrg/bp_updater.jl
new file mode 100644
index 00000000..ff247cc1
--- /dev/null
+++ b/src/beliefpropagationdmrg/bp_updater.jl
@@ -0,0 +1,46 @@
+using ITensors: contract
+using ITensors.ContractionSequenceOptimization: optimal_contraction_sequence
+using KrylovKit: eigsolve
+
+function default_krylov_kwargs()
+  return (; tol=1e-15, krylovdim=100, maxiter=100, verbosity=0, eager=false, ishermitian=true)
+end
+
+#TODO: Put inv_sqrt_mts onto ∂ψOψ_bpc_∂r beforehand. Need to do this in an efficient way without
+#precontracting ∂ψOψ_bpc_∂r and getting the index logic too messy
+function get_new_state(
+  ∂ψOψ_bpc_∂rs::Vector,
+  inv_sqrt_mts,
+  sqrt_mts,
+  state::ITensor;
+  sequences=["automatic" for i in length(∂ψOψ_bpc_∂rs)],
+)
+  state = noprime(contract([state; inv_sqrt_mts]))
+  states = ITensor[
+    noprime(contract([copy(state); ∂ψOψ_bpc_∂r]; sequence)) for
+    (∂ψOψ_bpc_∂r, sequence) in zip(∂ψOψ_bpc_∂rs, sequences)
+  ]
+  state = reduce(+, states)
+  return dag(noprime(contract([state; (inv_sqrt_mts)])))
+end
+
+function bp_eigsolve_updater(
+  init::ITensor,
+  ∂ψOψ_bpc_∂rs::Vector,
+  sqrt_mts,
+  inv_sqrt_mts;
+  krylov_kwargs=default_krylov_kwargs(),
+)
+  gauged_init = noprime(contract([copy(init); sqrt_mts]))
+  sequences = [
+    optimal_contraction_sequence([gauged_init; ∂ψOψ_bpc_∂r]) for ∂ψOψ_bpc_∂r in ∂ψOψ_bpc_∂rs
+  ]
+  get_new_state_ =
+    state -> get_new_state(∂ψOψ_bpc_∂rs, inv_sqrt_mts, sqrt_mts, state; sequences)
+  howmany = 1
+
+  vals, vecs, info = eigsolve(get_new_state_, gauged_init, howmany, :SR; krylov_kwargs...)
+  state = noprime(contract([first(vecs); inv_sqrt_mts]))
+
+  return state, first(vals)
+end
diff --git a/src/beliefpropagationdmrg/main.jl b/src/beliefpropagationdmrg/main.jl
new file mode 100644
index 00000000..8a7094c5
--- /dev/null
+++ b/src/beliefpropagationdmrg/main.jl
@@ -0,0 +1,34 @@
+using ITensorNetworks: random_tensornetwork
+using NamedGraphs.NamedGraphGenerators: named_comb_tree
+using ITensors: expect
+using Random
+
+include("graphsextensions.jl")
+include("tensornetworkoperators.jl")
+include("bp_dmrg.jl")
+include("bp_dmrg_V2.jl")
+include("utils.jl")
+
+Random.seed!(5634)
+
+function main()
+  g = named_grid((24, 1); periodic=true)
+  L = length(vertices(g))
+  h, hlongitudinal, J = 0.6, 0.2, 1.0
+  s = siteinds("S=1/2", g)
+  χ = 3
+  ψ0 = random_tensornetwork(s; link_space=χ)
+
+  H = ising(s; h, hl=hlongitudinal, J1=J)
+  #tnos = opsum_to_tnos_V2(s, H)
+  #tno = reduce(+, tnos)
+
+  energy_calc_fun =
+    (tn, bp_cache) -> sum(expect(tn, H; alg="bp", (cache!)=Ref(bp_cache))) / L
+
+  #ψfinal, energies = bp_dmrg(ψ0, tnos; no_sweeps=5, energy_calc_fun)
+  ψfinal, energies = bp_dmrg(ψ0, H; no_sweeps=5)
+  return final_mags = expect(ψfinal, "Z", ; alg="bp")
+end
+
+main()
diff --git a/src/caches/beliefpropagationcache.jl b/src/caches/beliefpropagationcache.jl
index b980a52f..ceb357df 100644
--- a/src/caches/beliefpropagationcache.jl
+++ b/src/caches/beliefpropagationcache.jl
@@ -14,13 +14,14 @@ using NamedGraphs.PartitionedGraphs:
 using SimpleTraits: SimpleTraits, Not, @traitfn
 
 default_message(inds_e) = ITensor[denseblocks(delta(i)) for i in inds_e]
+#default_message(inds_e) = ITensor[denseblocks(delta(inds_e))]  
 default_messages(ptn::PartitionedGraph) = Dictionary()
 default_message_norm(m::ITensor) = norm(m)
-function default_message_update(contract_list::Vector{ITensor}; kwargs...)
+function default_message_update(contract_list::Vector{ITensor}; normalize=true, kwargs...)
   sequence = optimal_contraction_sequence(contract_list)
   updated_messages = contract(contract_list; sequence, kwargs...)
   message_norm = norm(updated_messages)
-  if !iszero(message_norm)
+  if !iszero(message_norm) && normalize
     updated_messages /= message_norm
   end
   return ITensor[updated_messages]
@@ -132,6 +133,10 @@ function set_messages(cache::BeliefPropagationCache, messages)
   )
 end
 
+function delete_messages(cache::BeliefPropagationCache)
+  return BeliefPropagationCache(partitioned_tensornetwork(cache))
+end
+
 function environment(
   bp_cache::BeliefPropagationCache,
   partition_vertices::Vector{<:PartitionVertex};
@@ -157,6 +162,15 @@ function environment(bp_cache::BeliefPropagationCache, verts::Vector)
   return vcat(messages, central_tensors)
 end
 
+function factors(bp_cache::BeliefPropagationCache, vertices::Vector)
+  tn = tensornetwork(bp_cache)
+  return ITensor[tn[vertex] for vertex in vertices]
+end
+
+function factor(bp_cache::BeliefPropagationCache, vertex)
+  return only(factors(bp_cache, [vertex]))
+end
+
 function factor(bp_cache::BeliefPropagationCache, vertex::PartitionVertex)
   ptn = partitioned_tensornetwork(bp_cache)
   return collect(eachtensor(subgraph(ptn, vertex)))
@@ -173,6 +187,7 @@ function update_message(
 )
   vertex = src(edge)
   messages = environment(bp_cache, vertex; ignore_edges=PartitionEdge[reverse(edge)])
+
   state = factor(bp_cache, vertex)
 
   return message_update(ITensor[messages; state]; message_update_kwargs...)
@@ -258,6 +273,7 @@ Update the tensornetwork inside the cache
 """
 function update_factors(bp_cache::BeliefPropagationCache, factors)
   bp_cache = copy(bp_cache)
+  factors = copy(factors)
   tn = tensornetwork(bp_cache)
   for vertex in eachindex(factors)
     # TODO: Add a check that this preserves the graph structure.
@@ -309,3 +325,21 @@ end
 function scalar_factors_quotient(bp_cache::BeliefPropagationCache)
   return vertex_scalars(bp_cache), edge_scalars(bp_cache)
 end
+
+function ITensors.scalar(bp_cache::BeliefPropagationCache)
+  v_scalars, e_scalars = vertex_scalars(bp_cache), edge_scalars(bp_cache)
+  return prod(v_scalars) / prod(e_scalars)
+end
+
+#Renormalize all messages such that <m_{e}m_{reverse(e)}> = 1 on all edges
+function renormalize_messages(bp_cache::BeliefPropagationCache)
+  bp_cache = copy(bp_cache)
+  mts = messages(bp_cache)
+  for pe in partitionedges(partitioned_tensornetwork(bp_cache))
+    n = region_scalar(bp_cache, pe)
+    me, mer = only(mts[pe]), only(mts[reverse(pe)])
+    set!(mts, pe, ITensor[(1 / sqrt(n)) * me])
+    set!(mts, reverse(pe), ITensor[(1 / sqrt(n)) * mer])
+  end
+  return bp_cache
+end
diff --git a/src/expect.jl b/src/expect.jl
index e1d46a9f..a9fdfdaf 100644
--- a/src/expect.jl
+++ b/src/expect.jl
@@ -4,6 +4,12 @@ using ITensorMPS: ITensorMPS, expect
 
 default_expect_alg() = "bp"
 
+function ITensorMPS.expect(
+  ψ::AbstractITensorNetwork, ops; alg=default_expect_alg(), kwargs...
+)
+  return expect(Algorithm(alg), ψ, ops; kwargs...)
+end
+
 function ITensorMPS.expect(
   ψIψ::AbstractFormNetwork, op::Op; contract_kwargs=(; sequence="automatic"), kwargs...
 )
@@ -18,6 +24,27 @@ function ITensorMPS.expect(
   return numerator / denominator
 end
 
+function ITensorMPS.expect(
+  ψIψ::AbstractFormNetwork,
+  ops::Scaled{ComplexF64,Prod{Op}};
+  contract_kwargs=(; sequence="automatic"),
+  kwargs...,
+)
+  scalar = first(ops.args)
+  iszero(scalar) && return 0.0
+  n_ops = length(ops)
+  vs = site.(ops[1:n_ops])
+  op_strings = which_op.(ops[1:n_ops])
+  ψIψ_vs = [ψIψ[operator_vertex(ψIψ, v)] for v in vs]
+  sinds = [commonind(ψIψ[ket_vertex(ψIψ, v)], ψIψ_vs[i]) for (i, v) in enumerate(vs)]
+  operators = [ITensors.op(op_strings[i], sinds[i]) for i in 1:n_ops]
+  ∂ψIψ_∂v = environment(ψIψ, operator_vertices(ψIψ, vs); kwargs...)
+  numerator = contract(vcat(∂ψIψ_∂v, operators); contract_kwargs...)[]
+  denominator = contract(vcat(∂ψIψ_∂v, ψIψ_vs); contract_kwargs...)[]
+
+  return scalar * numerator / denominator
+end
+
 function ITensorMPS.expect(
   alg::Algorithm,
   ψ::AbstractITensorNetwork,
diff --git a/src/lib/ModelHamiltonians/src/ModelHamiltonians.jl b/src/lib/ModelHamiltonians/src/ModelHamiltonians.jl
index f54faea3..53bf6adf 100644
--- a/src/lib/ModelHamiltonians/src/ModelHamiltonians.jl
+++ b/src/lib/ModelHamiltonians/src/ModelHamiltonians.jl
@@ -108,8 +108,8 @@ heisenberg(N::Integer; kwargs...) = heisenberg(path_graph(N); kwargs...)
 """
 Next-to-nearest-neighbor Ising model (ZZX) on a general graph
 """
-function ising(g::AbstractGraph; J1=-1, J2=0, h=0)
-  (; J1, J2, h) = map(to_callable, (; J1, J2, h))
+function ising(g::AbstractGraph; J1=-1, J2=0, h=0, hl=0)
+  (; J1, J2, h, hl) = map(to_callable, (; J1, J2, h, hl))
   ℋ = OpSum()
   for e in edges(g)
     ℋ += J1(e), "Sz", src(e), "Sz", dst(e)
@@ -129,6 +129,7 @@ function ising(g::AbstractGraph; J1=-1, J2=0, h=0)
   end
   for v in vertices(g)
     ℋ += h(v), "Sx", v
+    ℋ += hl(v), "Sz", v
   end
   return ℋ
 end
diff --git a/src/solvers/insert/insert.jl b/src/solvers/insert/insert.jl
index a0250c95..9e0005c3 100644
--- a/src/solvers/insert/insert.jl
+++ b/src/solvers/insert/insert.jl
@@ -22,7 +22,7 @@ function default_inserter(
     indsTe = inds(state[ortho_vert])
     L, phi, spec = factorize(phi, indsTe; tags=tags(state, e), maxdim, mindim, cutoff)
     state[ortho_vert] = L
-
+    
   else
     v = ortho_vert
   end