diff --git a/Project.toml b/Project.toml
index 6f8f559c..7a5bc756 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "ITensorNetworks"
 uuid = "2919e153-833c-4bdc-8836-1ea460a35fc7"
-version = "0.21.1"
+version = "0.22.0-DEV"
 authors = ["Matthew Fishman <mfishman@flatironinstitute.org>, Joseph Tindall <jtindall@flatironinstitute.org> and contributors"]
 
 [workspace]
@@ -15,7 +15,6 @@ Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"
 ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
 DataGraphs = "b5a273c3-7e6c-41f6-98bd-8d7f1525a36a"
 Dictionaries = "85a47980-9c8c-11e8-2b9f-f7ca1fa99fb4"
-Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
 ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5"
@@ -54,7 +53,6 @@ Compat = "3, 4"
 ConstructionBase = "1.6"
 DataGraphs = "0.4"
 Dictionaries = "0.4"
-Distributions = "0.25.86"
 DocStringExtensions = "0.9"
 Graphs = "1.8"
 ITensors = "0.7, 0.8, 0.9"
diff --git a/docs/Project.toml b/docs/Project.toml
index 35064770..bf58fae3 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -16,7 +16,7 @@ path = ".."
 Documenter = "1.10"
 Graphs = "1"
 ITensorFormatter = "0.2.27"
-ITensorNetworks = "0.21"
+ITensorNetworks = "0.22"
 ITensors = "0.9"
 Literate = "2.20.1"
 NamedGraphs = "0.11"
diff --git a/docs/make.jl b/docs/make.jl
index 18cf0162..0e4e7e01 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -14,7 +14,7 @@ DocMeta.setdocmeta!(
         using Graphs: dst, edges, src, vertices
         using ITensorNetworks
         using ITensorNetworks:
-            TreeTensorNetwork, expect, loginner, mps, orthogonalize, siteinds, truncate, ttn
+            TreeTensorNetwork, expect, loginner, orthogonalize, siteinds, truncate
         using ITensors: inner
         using LinearAlgebra: norm, normalize
         using OMEinsumContractionOrders
diff --git a/docs/src/computing_properties.md b/docs/src/computing_properties.md
index f2b221ed..71396880 100644
--- a/docs/src/computing_properties.md
+++ b/docs/src/computing_properties.md
@@ -2,7 +2,8 @@
 
 ```@setup main
 using Graphs: edges, vertices
-using ITensorNetworks: ITensorNetwork, expect, inner, loginner, normalize, siteinds, ttn
+using ITensorNetworks:
+    ITensorNetwork, TreeTensorNetwork, expect, inner, loginner, normalize, siteinds
 using ITensors: Index, random_itensor
 using LinearAlgebra: norm
 using NamedGraphs.GraphsExtensions: incident_edges
@@ -20,10 +21,10 @@ end
 
 g = named_grid((4,))
 s = siteinds("S=1/2", g)
-phi = normalize(ttn(random_state(g, s; link_space = 2)))
-psi = normalize(ttn(random_state(g, s; link_space = 2)))
-x = normalize(ttn(random_state(g, s; link_space = 2)))
-y = normalize(ttn(random_state(g, s; link_space = 2)))
+phi = normalize(TreeTensorNetwork(random_state(g, s; link_space = 2)))
+psi = normalize(TreeTensorNetwork(random_state(g, s; link_space = 2)))
+x = normalize(TreeTensorNetwork(random_state(g, s; link_space = 2)))
+y = normalize(TreeTensorNetwork(random_state(g, s; link_space = 2)))
 v = first(vertices(psi))
 ```
 
diff --git a/docs/src/deprecated_methods.md b/docs/src/deprecated_methods.md
index 48dd64bc..35a0b5dc 100644
--- a/docs/src/deprecated_methods.md
+++ b/docs/src/deprecated_methods.md
@@ -2,43 +2,6 @@
 
 Suggestions of methods which could be deleted.
 
-## ITensorNetwork Methods
-
-#### ITensorNetwork Constructors
-
-* From a named graph, forwards to construction from `IndsNetwork` (`itensornetwork.jl`):
-  ```julia
-  ITensorNetwork{V}(g::NamedGraph)
-  ITensorNetwork(eltype::Type, undef::UndefInitializer, graph::AbstractNamedGraph; kws...)
-  ITensorNetwork(f, graph::AbstractNamedGraph; kwargs...)
-  ITensorNetwork(graph::AbstractNamedGraph; kwargs...)
-  ```
-
-* From a simple graph, forwards to construction from `IndsNetwork` (`itensornetwork.jl`):
-  ```julia
-  ITensorNetwork(eltype::Type, undef::UndefInitializer, graph::AbstractSimpleGraph; kws...)
-  ITensorNetwork(f, graph::AbstractSimpleGraph; kwargs...)
-  ITensorNetwork(graph::AbstractSimpleGraph; kwargs...)
-  ```
-
-* From a function over vertices or from a "value" (e.g. a string like `"Up"`,
-  an `Op`, an array, or a per-vertex dict/array) that is converted to a callable and used
-  to initialize each vertex tensor (`itensornetwork.jl`):
-  ```julia
-  ITensorNetwork(value, is::IndsNetwork; kwargs...)
-  ITensorNetwork(elt::Type, f, is::IndsNetwork; link_space = trivial_space(is), kws...)
-  ITensorNetwork(itensor_constructor::Function, is::IndsNetwork; link_space = trivial_space(is), kwargs...)
-  ```
-
-* Construct an `ITensorNetwork` from an `IndsNetwork`. Initializes ITensors with `undef` storage on each vertex
-  of the `IndsNetwork` with the corresponding indices (`itensornetwork.jl`):
-  ```julia
-  ITensorNetwork(eltype::Type, undef::UndefInitializer, is::IndsNetwork; kwargs...)
-  ITensorNetwork(eltype::Type, is::IndsNetwork; kwargs...)
-  ITensorNetwork(undef::UndefInitializer, is::IndsNetwork; kwargs...)
-  ITensorNetwork(is::IndsNetwork; kwargs...)
-  ```
-
 ## Global Operations on ITensorNetworks
 
 
diff --git a/docs/src/experimental_methods.md b/docs/src/experimental_methods.md
index 14d28373..1c1895fc 100644
--- a/docs/src/experimental_methods.md
+++ b/docs/src/experimental_methods.md
@@ -50,7 +50,7 @@ Methods which still need to be discussed, modified, or deprecated.
 
 * From an `OpSum`, using `opsum_to_ttn.jl` code:
   ```julia
-  ttn(os::OpSum, sites::IndsNetwork; kws...)
+  TreeTensorNetwork(os::OpSum, sites::IndsNetwork; kws...)
   ```
 
 * From `OpSum`, assuming path graph (`opsum_to_ttn.jl`):
diff --git a/docs/src/interface_methods.md b/docs/src/interface_methods.md
index 2c40c497..7eb40df5 100644
--- a/docs/src/interface_methods.md
+++ b/docs/src/interface_methods.md
@@ -6,20 +6,20 @@ Recommended methods for building applications on top of ITensorNetworks.
 
 These ITensorNetwork constructor interfaces are foundational to other constructors:
 
-* From dictionary-like objects, including other `ITensorNetwork` objects, or a `Dict`
-  from vertices to `ITensor`. The `tensors` can be a collection of `ITensor`s in which case
-  the vertex labels are auto-assigned to `eachindex(tensors)` and edges are inferred
-  from shared indices between the `ITensor`s.
+* From a collection of `ITensor`s indexed by vertex (a `Dict`, `Dictionary`, or
+  `Vector{ITensor}` with linear-index vertex labels). Edges are inferred from shared
+  `Index` objects between the tensors.
   ```julia
   # `keys(tensors)` are vertices, `values(tensors)` are tensors on those vertices
   ITensorNetwork(tensors)
   ITensorNetwork{V}(tensors)
   ```
 
-* From collections of vertices and tensors. E.g. a `Vector{Int}` and a `Vector{ITensor}`.
+* From a collection of `ITensor`s placed at the vertices of a given `NamedGraph`. No
+  edge inference; the graph's edges are used as-is.
   ```julia
-  ITensorNetwork(vertices, tensors)
-  ITensorNetwork{V}(vertices, tensors)
+  ITensorNetwork(tensors, graph::NamedGraph)
+  ITensorNetwork{V}(tensors, graph::NamedGraph)
   ```
 
 
diff --git a/docs/src/itensor_networks.md b/docs/src/itensor_networks.md
index 7b225c1f..d14dc7fc 100644
--- a/docs/src/itensor_networks.md
+++ b/docs/src/itensor_networks.md
@@ -11,17 +11,13 @@ Key facts:
 - The underlying graph is a [`NamedGraph`](https://github.com/ITensor/NamedGraphs.jl), so
   vertices can be any hashable Julia value: integers, tuples, strings, etc.
 - Each vertex holds exactly one `ITensor`.
-- Edges and link indices are either inferred from shared `Index` objects (when constructing
-  from a collection of `ITensor`s) or inserted automatically (when constructing from an
-  `IndsNetwork`).
+- Edges are either inferred from shared `Index` objects (when constructing from a
+  collection of `ITensor`s) or taken from a graph passed explicitly alongside the tensors.
 
 ## Construction
 
-The most common entry point is an `IndsNetwork` — a graph whose vertices and edges carry
-`Index` objects.  Generate site indices with the `siteinds` function which takes a site
-type string (such as "S=1/2" or "Electron") and a NamedGraph. The NamedGraph can be 
-generated from functions such as `named_grid`, `named_comb_tree`, etc. from the NamedGraphs.jl
-`NamedGraphGenerators` module:
+When you already have `ITensor`s in hand, edges are inferred automatically from shared
+indices:
 
 ```@example main
 using Graphs: edges, ne, neighbors, nv, vertices
@@ -29,31 +25,44 @@ using ITensorNetworks: ITensorNetwork, add, linkinds, siteinds
 using ITensors: Index, ITensor
 using NamedGraphs.NamedGraphGenerators: named_grid
 
-# 3×3 square-lattice tensor network
-g = named_grid((3, 3))
-s = siteinds("S=1/2", g)  # one spin-½ Index per vertex
-
-# Zero-initialized, bond dimension 2
-ψ = ITensorNetwork(s; link_space = 2)
-
-# Product state — every site in the |↑⟩ state
-ψ = ITensorNetwork("Up", s)
+i, j, k = Index(2, "i"), Index(2, "j"), Index(2, "k")
+A, B, C = ITensor(i, j), ITensor(j, k), ITensor(k)
 
-# Staggered initialization with a vertex-dependent function
-ψ = ITensorNetwork(v -> isodd(sum(v)) ? "Up" : "Dn", s)
+tn = ITensorNetwork([A, B, C])  # integer vertices 1, 2, 3
+tn = ITensorNetwork(Dict("A" => A, "B" => B, "C" => C))  # named vertices via a Dict
 ```
 
-When you already have `ITensor`s in hand, edges are inferred automatically from shared
-indices:
+If you want to control edges directly — for example to build a structured network on a
+prescribed lattice and fill in tensors later — pass a `NamedGraph` along with a
+collection of `ITensor`s indexed by vertex. To create a tensor network with shared link
+indices on each edge, build the indices once per edge and reuse them at both endpoints:
 
 ```@example main
-i, j, k = Index(2, "i"), Index(2, "j"), Index(2, "k")
-A, B, C = ITensor(i, j), ITensor(j, k), ITensor(k)
+using ITensors: random_itensor
+using NamedGraphs: NamedGraph
+using NamedGraphs.GraphsExtensions: edgetype, incident_edges
 
-tn = ITensorNetwork([A, B, C])  # integer vertices 1, 2, 3
-tn = ITensorNetwork(Dict("A" => A, "B" => B, "C" => C))  # named vertices via a Dict
+g = NamedGraph(named_grid((3, 3)))
+s = siteinds("S=1/2", g)  # one spin-½ site Index per vertex
+
+# One shared link Index per edge; bond dimension χ
+χ = 2
+links = Dict(e => Index(χ, "Link") for e in edges(g))
+
+# Per-vertex tensor: the site Index plus the link Index of every incident edge
+tensors = Dict(map(collect(vertices(g))) do v
+    site_v = s[v]
+    link_v = [haskey(links, e) ? links[e] : links[reverse(e)] for e in incident_edges(g, v)]
+    return v => random_itensor(site_v..., link_v...)
+end)
+
+ψ = ITensorNetwork(tensors, g)
 ```
 
+Higher-level construction routines (random networks, product states, OpSum-derived
+TTNs, etc.) are provided by sibling functions like `TreeTensorNetwork(opsum, sites)`
+and the test-only helpers in `test/utils.jl`.
+
 ```@docs; canonical=false
 ITensorNetworks.ITensorNetwork
 ```
@@ -77,7 +86,7 @@ linkinds(ψ)  # IndsNetwork of bond (virtual) indices
 Two networks with the same graph and site indices can be added. The result represents the
 tensor network `ψ₁ + ψ₂` and has bond dimension equal to the **sum** of the two input bond
 dimensions. Individual bonds of the result can be recompressed with `truncate(tn, edge)`.
-For `TreeTensorNetwork`, the no-argument form `truncate(ttn; kwargs...)` sweeps and
+For `TreeTensorNetwork`, the no-argument form `truncate(tn; kwargs...)` sweeps and
 recompresses all bonds at once.
 
 ```@example main
@@ -101,7 +110,7 @@ edge = (1, 2) => (1, 3)
 ```
 
 Truncation parameters (`cutoff`, `maxdim`, `mindim`, …) are forwarded to `ITensors.svd`.
-For a `TreeTensorNetwork`, the sweep-based `truncate(ttn; kwargs...)` is usually more
+For a `TreeTensorNetwork`, the sweep-based `truncate(tn; kwargs...)` is usually more
 convenient because it recompresses the entire network at once with controlled errors;
 see the [Tree Tensor Networks](@ref) page.
 
diff --git a/docs/src/solvers.md b/docs/src/solvers.md
index 34f50c2c..d8bf4976 100644
--- a/docs/src/solvers.md
+++ b/docs/src/solvers.md
@@ -17,7 +17,8 @@ variational sweep algorithm.
 
 ```@example main
 using Graphs: vertices
-using ITensorNetworks: ITensorNetwork, dmrg, dst, edges, normalize, siteinds, src, ttn
+using ITensorNetworks:
+    ITensorNetwork, TreeTensorNetwork, dmrg, dst, edges, normalize, siteinds, src
 using ITensors: Index, OpSum, random_itensor
 using NamedGraphs.GraphsExtensions: incident_edges
 using NamedGraphs.NamedGraphGenerators: named_comb_tree
@@ -29,7 +30,7 @@ function random_state(g, s; link_space)
         v => random_itensor(only(s[v]), (l[e] for e in incident_edges(g, v))...)
             for v in vertices(g)
     )
-    return ITensorNetwork(ts)
+    return ITensorNetwork(ts, g)
 end
 
 # Build a Heisenberg Hamiltonian on a comb tree
@@ -41,11 +42,11 @@ H = let h = OpSum()
         h += 0.5, "S-", src(e), "S+", dst(e)
         h += "Sz", src(e), "Sz", dst(e)
     end
-    ttn(h, s)
+    TreeTensorNetwork(h, s)
 end
 
 # Random initial state (normalise first!)
-psi0 = normalize(ttn(random_state(g, s; link_space = 2)))
+psi0 = normalize(TreeTensorNetwork(random_state(g, s; link_space = 2)))
 
 # Run DMRG
 energy, psi = dmrg(H, psi0;
diff --git a/docs/src/tree_tensor_networks.md b/docs/src/tree_tensor_networks.md
index 52ed2231..ff740ffb 100644
--- a/docs/src/tree_tensor_networks.md
+++ b/docs/src/tree_tensor_networks.md
@@ -11,39 +11,20 @@ records which vertices currently form the orthogonality center of the network. A
 update this field as the gauge changes.
 
 **MPS** (matrix product states) are the special case of a `TreeTensorNetwork` on a
-1D path graph. The [`mps`](@ref ITensorNetworks.mps) constructor enforces this topology
-and provides a convenient interface for 1D calculations.
+1D path graph.
 
 ## Construction
 
-### From an `IndsNetwork` or graph
+### From an `OpSum` (Hamiltonian)
 
-```@example main
-using Graphs: vertices
-using ITensorNetworks: ITensorNetwork, TreeTensorNetwork, mps, ortho_region, orthogonalize,
-    siteinds, ttn
-using ITensors: ITensors
-using LinearAlgebra: norm
-using NamedGraphs.NamedGraphGenerators: named_comb_tree
-
-# Comb-tree TTN (a popular tree topology for 2D-like systems)
-g = named_comb_tree((3, 2))
-sites = siteinds("S=1/2", g)
-
-psi = ttn(sites)  # zero-initialised
-psi = ttn(v -> "Up", sites)  # product state
-
-# 1D MPS
-s1d = siteinds("S=1/2", 6)
-mps_state = mps(v -> "Up", s1d)  # product MPS
-```
+A common way to obtain a Hamiltonian-shaped TTN is to convert an `OpSum` over an
+`IndsNetwork` of site indices.
 
 ```@docs; canonical=false
-ITensorNetworks.ttn
-ITensorNetworks.mps
+ITensorNetworks.TreeTensorNetwork(::ITensors.Ops.OpSum, ::ITensorNetworks.IndsNetwork)
 ```
 
-### The `TreeTensorNetwork` type and conversion
+### From an existing `ITensorNetwork`
 
 The `TreeTensorNetwork` struct wraps an `ITensorNetwork` and records the current
 orthogonality region. Use the `TreeTensorNetwork` constructor to convert a plain
@@ -51,8 +32,35 @@ orthogonality region. Use the `TreeTensorNetwork` constructor to convert a plain
 gauge metadata when you need a plain network again.
 
 ```@example main
-itn = ITensorNetwork(psi)  # TTN → ITensorNetwork
-psi = TreeTensorNetwork(itn)  # ITensorNetwork → TTN
+using Graphs: edges, vertices
+using ITensorNetworks: ITensorNetwork, TreeTensorNetwork, ortho_region, orthogonalize,
+    siteinds
+using ITensors: ITensors, Index, random_itensor
+using LinearAlgebra: norm
+using NamedGraphs: NamedGraph
+using NamedGraphs.GraphsExtensions: incident_edges
+using NamedGraphs.NamedGraphGenerators: named_comb_tree
+
+# Comb-tree TTN (a popular tree topology for 2D-like systems)
+g = NamedGraph(named_comb_tree((3, 2)))
+sites = siteinds("S=1/2", g)
+
+# Build a structured `ITensorNetwork` with shared link indices on each edge
+χ = 2
+links = Dict(e => Index(χ, "Link") for e in edges(g))
+tensors = Dict(map(collect(vertices(g))) do v
+    site_v = sites[v]
+    link_v = [haskey(links, e) ? links[e] : links[reverse(e)] for e in incident_edges(g, v)]
+    return v => random_itensor(site_v..., link_v...)
+end)
+itn = ITensorNetwork(tensors, g)
+psi = TreeTensorNetwork(itn)
+```
+
+To strip the gauge metadata back to a plain `ITensorNetwork`:
+
+```@example main
+itn_again = ITensorNetwork(psi)  # TTN → ITensorNetwork
 ```
 
 ```@docs; canonical=false
@@ -67,14 +75,16 @@ tree edges. Truncation parameters (e.g. `cutoff`, `maxdim`) are forwarded to the
 factorisation step.
 
 ```@example main
-g = named_comb_tree((3, 1))
-sites = siteinds("S=1/2", g)
-A = ITensors.random_itensor(only(sites[(1, 1)]), only(sites[(2, 1)]), only(sites[(3, 1)]))
-ttn_A = ttn(A, sites)
+g_small = named_comb_tree((3, 1))
+sites_small = siteinds("S=1/2", g_small)
+A = ITensors.random_itensor(
+    only(sites_small[(1, 1)]), only(sites_small[(2, 1)]), only(sites_small[(3, 1)])
+)
+ttn_A = TreeTensorNetwork(A, sites_small)
 ```
 
 ```@docs; canonical=false
-ITensorNetworks.ttn(::ITensors.ITensor, ::ITensorNetworks.IndsNetwork)
+ITensorNetworks.TreeTensorNetwork(::ITensors.ITensor, ::ITensorNetworks.IndsNetwork)
 ```
 
 ## Orthogonal Gauge
diff --git a/examples/Project.toml b/examples/Project.toml
index 3a9411e4..7899d300 100644
--- a/examples/Project.toml
+++ b/examples/Project.toml
@@ -5,4 +5,4 @@ ITensorNetworks = "2919e153-833c-4bdc-8836-1ea460a35fc7"
 path = ".."
 
 [compat]
-ITensorNetworks = "0.21"
+ITensorNetworks = "0.22"
diff --git a/src/abstractindsnetwork.jl b/src/abstractindsnetwork.jl
index a5799703..49af846e 100644
--- a/src/abstractindsnetwork.jl
+++ b/src/abstractindsnetwork.jl
@@ -1,4 +1,4 @@
-using .ITensorsExtensions: ITensorsExtensions, promote_indtype
+using .ITensorsExtensions: ITensorsExtensions, promote_indtype, trivial_space
 using DataGraphs: DataGraphs, AbstractDataGraph, DataGraph, IsUnderlyingGraph, edge_data,
     get_edge_data, get_vertex_data, is_edge_assigned, is_vertex_assigned, map_data,
     set_edge_data!, set_vertex_data!, underlying_graph_type, vertex_data
diff --git a/src/abstractitensornetwork.jl b/src/abstractitensornetwork.jl
index 11f91c96..83eacb51 100644
--- a/src/abstractitensornetwork.jl
+++ b/src/abstractitensornetwork.jl
@@ -1,14 +1,14 @@
-using .ITensorsExtensions: ITensorsExtensions, indtype, promote_indtype
+using .ITensorsExtensions: ITensorsExtensions, indtype, promote_indtype, trivial_space
 using Adapt: Adapt, adapt, adapt_structure
 using DataGraphs: DataGraphs, edge_data, get_vertex_data, is_vertex_assigned,
     set_vertex_data!, underlying_graph, underlying_graph_type, vertex_data
 using Dictionaries: Dictionary
 using Graphs: Graphs, Graph, add_edge!, add_vertex!, bfs_tree, center, dst, edges, edgetype,
-    ne, neighbors, rem_edge!, src, vertices
+    has_edge, ne, neighbors, rem_edge!, src, vertices
 using ITensors: ITensors, @Algorithm_str, ITensor, addtags, combiner, commoninds,
     commontags, contract, dag, inds, noprime, onehot, prime, replaceprime, replacetags,
     setprime, settags, sim, swaptags, tags
-using LinearAlgebra: LinearAlgebra, factorize
+using LinearAlgebra: LinearAlgebra, factorize, qr
 using MacroTools: @capture
 using NDTensors: NDTensors, Algorithm, dim, scalartype
 using NamedGraphs.GraphsExtensions:
@@ -136,7 +136,7 @@ function fix_edges!(tn::AbstractITensorNetwork, v)
         if v ≠ vertex
             edge = v => vertex
             if !isempty(linkinds(tn, edge))
-                add_edge!(tn, edge)
+                add_edge!(data_graph(tn), edge)
             end
         end
     end
@@ -181,7 +181,7 @@ function Base.union(tn1::AbstractITensorNetwork, tn2::AbstractITensorNetwork; kw
     for v1 in vertices(tn1)
         for v2 in vertices(tn2)
             if !isempty(linkinds(tn, v1 => v2))
-                add_edge!(tn, v1 => v2)
+                add_edge!(data_graph(tn), v1 => v2)
             end
         end
     end
@@ -274,7 +274,7 @@ end
 #
 
 function _siteinds(tn::AbstractITensorNetwork, vertex)
-    s = inds(tn[vertex])
+    s = collect(inds(tn[vertex]))
     for v in neighbors(tn, vertex)
         s = setdiff(s, inds(tn[v]))
     end
@@ -471,10 +471,10 @@ function NDTensors.contract(
     rem_vertex!(tn, dst(edge))
     add_vertex!(tn, merged_vertex)
     for n_src in neighbors_src
-        add_edge!(tn, merged_vertex => n_src)
+        add_edge!(data_graph(tn), merged_vertex => n_src)
     end
     for n_dst in neighbors_dst
-        add_edge!(tn, merged_vertex => n_dst)
+        add_edge!(data_graph(tn), merged_vertex => n_dst)
     end
     @preserve_graph tn[merged_vertex] = new_itensor
     return tn
@@ -561,11 +561,11 @@ function LinearAlgebra.factorize(
     add_vertex!(tn, X_vertex)
     add_vertex!(tn, Y_vertex)
 
-    add_edge!(tn, X_vertex => Y_vertex)
+    add_edge!(data_graph(tn), X_vertex => Y_vertex)
     for nX in neighbors_X
-        add_edge!(tn, X_vertex => nX)
+        add_edge!(data_graph(tn), X_vertex => nX)
     end
-    add_edge!(tn, Y_vertex => dst(edge))
+    add_edge!(data_graph(tn), Y_vertex => dst(edge))
     @preserve_graph tn[X_vertex] = X
     @preserve_graph tn[Y_vertex] = Y
     return tn
@@ -681,7 +681,7 @@ Truncation parameters are passed as keyword arguments and forwarded to `ITensors
   - `mindim`: Minimum number of singular values to keep.
 
 This operates on a single bond. For `TreeTensorNetwork`, the no-argument form
-`truncate(ttn; kwargs...)` sweeps all bonds and is generally preferred for full
+`truncate(tn; kwargs...)` sweeps all bonds and is generally preferred for full
 recompression after addition or subspace expansion.
 
 See also: `Base.truncate(::AbstractTreeTensorNetwork)`.
@@ -855,6 +855,27 @@ function insert_linkinds(
     return tn
 end
 
+# Factor `tn[src(edge)]` via `f` along the inds it doesn't share with
+# `tn[dst(edge)]`, leaving the new factor on `src(edge)` and contracting
+# the residual into `tn[dst(edge)]`. With `f = qr` and no shared inds yet,
+# this threads a fresh link Index between the two endpoints. With shared
+# inds present, it regauges across the existing link.
+function factorize_edge!(f, tn::AbstractITensorNetwork, edge)
+    src_t, dst_t = tn[src(edge)], tn[dst(edge)]
+    left_inds = setdiff(inds(src_t), inds(dst_t))
+    x, y = f(src_t, left_inds)
+    @preserve_graph tn[src(edge)] = x
+    @preserve_graph tn[dst(edge)] = y * dst_t
+    return tn
+end
+
+function Graphs.add_edge!(tn::AbstractITensorNetwork, edge)
+    has_edge(tn, edge) && return false
+    add_edge!(data_graph(tn), edge)
+    factorize_edge!(qr, tn, edge)
+    return true
+end
+
 # TODO: What to output? Could be an `IndsNetwork`. Or maybe
 # that would be a different function `commonindsnetwork`.
 # Even in that case, this could output a `Dictionary`
diff --git a/src/formnetworks/bilinearformnetwork.jl b/src/formnetworks/bilinearformnetwork.jl
index e0bc24aa..4c632c1a 100644
--- a/src/formnetworks/bilinearformnetwork.jl
+++ b/src/formnetworks/bilinearformnetwork.jl
@@ -84,14 +84,16 @@ function BilinearFormNetwork(
     bra_site_inds = mapreduce(v -> siteinds(bra, v), vcat, vertices(bra); init = Index[])
     ket_site_inds = mapreduce(v -> siteinds(ket, v), vcat, vertices(ket); init = Index[])
     @assert issetequal(bra_site_inds, ket_site_inds)
-    link_space = isempty(bra_site_inds) ? 1 : nothing
     s = siteinds(ket)
     s_mapped = dual_site_index_map(s)
     operator_inds = union_all_inds(s, s_mapped)
 
-    O = ITensorNetwork(operator_inds; link_space) do v
-        return inds -> itensor_identity_map(scalartype(ket), s[v] .=> s_mapped[v])
+    g = NamedGraph(underlying_graph(operator_inds))
+    ts = Dict{vertextype(g), ITensor}()
+    for v in vertices(operator_inds)
+        ts[v] = itensor_identity_map(scalartype(ket), s[v] .=> s_mapped[v])
     end
+    O = ITensorNetwork(ts, g)
     O = adapt(promote_type(datatype(bra), datatype(ket)), O)
     return BilinearFormNetwork(O, bra, ket; dual_site_index_map, kwargs...)
 end
diff --git a/src/indsnetwork.jl b/src/indsnetwork.jl
index c7aadb3e..c7f6c851 100644
--- a/src/indsnetwork.jl
+++ b/src/indsnetwork.jl
@@ -302,15 +302,3 @@ Base.copy(is::IndsNetwork) = IndsNetwork(copy(data_graph(is)))
 function map_inds(f, is::IndsNetwork, args...; sites = nothing, links = nothing, kwargs...)
     return map_data(i -> f(i, args...; kwargs...), is; vertices = sites, edges = links)
 end
-
-#
-# Visualization
-#
-
-# TODO: Move to an `ITensorNetworksVisualizationInterfaceExt`
-# package extension (and define a `VisualizationInterface` package
-# based on `ITensorVisualizationCore`.).
-using ITensors.ITensorVisualizationCore: ITensorVisualizationCore, visualize
-function ITensorVisualizationCore.visualize(is::IndsNetwork, args...; kwargs...)
-    return visualize(ITensorNetwork(is), args...; kwargs...)
-end
diff --git a/src/itensornetwork.jl b/src/itensornetwork.jl
index 7073c5da..e0f6ce38 100644
--- a/src/itensornetwork.jl
+++ b/src/itensornetwork.jl
@@ -1,7 +1,5 @@
-using .ITensorsExtensions: trivial_space
 using DataGraphs: DataGraphs, DataGraph
-using Dictionaries: Indices
-using ITensors: ITensors, ITensor, op
+using ITensors: ITensors, ITensor
 using NamedGraphs: NamedGraphs, NamedEdge, NamedGraph, similar_graph, vertextype
 
 """
@@ -29,42 +27,18 @@ ITensorNetwork(tensors)
 ITensorNetwork(tensors, graph::NamedGraph)
 ```
 
-**From an `IndsNetwork`** (most common in legacy code):
-
-```julia
-ITensorNetwork(is::IndsNetwork; link_space = 1)
-ITensorNetwork(f, is::IndsNetwork; link_space = 1)
-ITensorNetwork(eltype, undef, is::IndsNetwork; link_space = 1)
-```
-
-  - With no function argument `f`, tensors are initialized to zero.
-  - With a function `f(v)` that returns a state label (e.g. `"Up"`, `"Dn"`) or
-    an `ITensor` constructor, tensors are initialized accordingly.
-  - `link_space` controls the bond-index dimension (default 1).
-
-**From a graph (site indices inferred as trivial):**
-
-```julia
-ITensorNetwork(graph::AbstractNamedGraph; link_space = 1)
-ITensorNetwork(f, graph::AbstractNamedGraph; link_space = 1)
-```
-
 # Example
 
 ```jldoctest
-julia> using NamedGraphs.NamedGraphGenerators: named_grid
+julia> using ITensors: Index, ITensor
 
-julia> g = named_grid((4,));
+julia> i, j, k = Index(2, "i"), Index(2, "j"), Index(2, "k");
 
-julia> s = siteinds("S=1/2", g);
-
-julia> tn = ITensorNetwork(s; link_space = 2);
-
-julia> tn = ITensorNetwork("Up", s);
+julia> tn = ITensorNetwork([ITensor(i, j), ITensor(j, k)]);
 
 ```
 
-See also: `IndsNetwork`, [`ttn`](@ref ITensorNetworks.ttn), [`TreeTensorNetwork`](@ref ITensorNetworks.TreeTensorNetwork).
+See also: `IndsNetwork`, [`TreeTensorNetwork`](@ref ITensorNetworks.TreeTensorNetwork).
 """
 const _ITensorCollection = Union{
     AbstractVector{<:ITensor},
@@ -147,9 +121,6 @@ function ITensorNetwork{V}(tn::ITensorNetwork) where {V}
     tensors = Dict{V, ITensor}(v => tn[v] for v in vertices(tn))
     return ITensorNetwork(tensors, g)
 end
-function ITensorNetwork{V}(g::NamedGraph) where {V}
-    return ITensorNetwork(NamedGraph{V}(g))
-end
 
 ITensorNetwork(tn::ITensorNetwork) = copy(tn)
 
@@ -167,153 +138,3 @@ function NamedGraphs.similar_graph(tn::ITensorNetwork, underlying_graph::Abstrac
     default = Dict{vertextype(g), ITensor}(v => ITensor() for v in vertices(g))
     return ITensorNetwork(default, g)
 end
-
-#
-# Construction from underyling named graph
-#
-
-function ITensorNetwork(
-        eltype::Type, undef::UndefInitializer, graph::AbstractNamedGraph; kwargs...
-    )
-    return ITensorNetwork(eltype, undef, IndsNetwork(graph; kwargs...))
-end
-
-function ITensorNetwork(f, graph::AbstractNamedGraph; kwargs...)
-    return ITensorNetwork(f, IndsNetwork(graph; kwargs...))
-end
-
-function ITensorNetwork(graph::AbstractNamedGraph; kwargs...)
-    return ITensorNetwork(IndsNetwork(graph; kwargs...))
-end
-
-#
-# Construction from underyling simple graph
-#
-
-function ITensorNetwork(
-        eltype::Type, undef::UndefInitializer, graph::AbstractSimpleGraph; kwargs...
-    )
-    return ITensorNetwork(eltype, undef, IndsNetwork(graph; kwargs...))
-end
-
-function ITensorNetwork(f, graph::AbstractSimpleGraph; kwargs...)
-    return ITensorNetwork(f, IndsNetwork(graph); kwargs...)
-end
-
-function ITensorNetwork(graph::AbstractSimpleGraph; kwargs...)
-    return ITensorNetwork(IndsNetwork(graph); kwargs...)
-end
-
-#
-# Construction from IndsNetwork
-#
-
-function ITensorNetwork(eltype::Type, undef::UndefInitializer, is::IndsNetwork; kwargs...)
-    return ITensorNetwork(is; kwargs...) do v
-        return (inds...) -> ITensor(eltype, undef, inds...)
-    end
-end
-
-function ITensorNetwork(eltype::Type, is::IndsNetwork; kwargs...)
-    return ITensorNetwork(is; kwargs...) do v
-        return (inds...) -> ITensor(eltype, inds...)
-    end
-end
-
-function ITensorNetwork(undef::UndefInitializer, is::IndsNetwork; kwargs...)
-    return ITensorNetwork(is; kwargs...) do v
-        return (inds...) -> ITensor(undef, inds...)
-    end
-end
-
-function ITensorNetwork(is::IndsNetwork; kwargs...)
-    return ITensorNetwork(is; kwargs...) do v
-        return (inds...) -> ITensor(inds...)
-    end
-end
-
-# TODO: Handle `eltype` and `undef` through `generic_state`.
-# `inds` are stored in a `NamedTuple`
-function generic_state(f, inds::NamedTuple)
-    return generic_state(f, reduce(vcat, inds.linkinds; init = inds.siteinds))
-end
-
-function generic_state(f, inds::Vector)
-    return f(inds)
-end
-function generic_state(a::AbstractArray, inds::Vector)
-    return itensor(a, inds)
-end
-function generic_state(x::Op, inds::NamedTuple)
-    # TODO: Figure out what to do if there is more than one site.
-    if !isempty(inds.siteinds)
-        @assert length(inds.siteinds) == 2
-        i = inds.siteinds[findfirst(i -> plev(i) == 0, inds.siteinds)]
-        @assert i' ∈ inds.siteinds
-        site_tensors = [op(x.which_op, i)]
-    else
-        site_tensors = []
-    end
-    link_tensors = [[onehot(i => 1) for i in inds.linkinds[e]] for e in keys(inds.linkinds)]
-    return contract(reduce(vcat, link_tensors; init = site_tensors))
-end
-function generic_state(s::AbstractString, inds::NamedTuple)
-    # TODO: Figure out what to do if there is more than one site.
-    site_tensors = [ITensors.state(s, only(inds.siteinds))]
-    link_tensors = [[onehot(i => 1) for i in inds.linkinds[e]] for e in keys(inds.linkinds)]
-    return contract(reduce(vcat, link_tensors; init = site_tensors))
-end
-
-# TODO: This is similar to `ModelHamiltonians.to_callable`,
-# try merging the two.
-to_callable(value::Type) = value
-to_callable(value::Function) = value
-function to_callable(value::AbstractDict)
-    return Base.Fix1(getindex, value) ∘ keytype(value)
-end
-function to_callable(value::AbstractDictionary)
-    return Base.Fix1(getindex, value) ∘ keytype(value)
-end
-function to_callable(value::AbstractArray{<:Any, N}) where {N}
-    return Base.Fix1(getindex, value) ∘ CartesianIndex
-end
-to_callable(value) = Returns(value)
-
-function ITensorNetwork(value, is::IndsNetwork; kwargs...)
-    return ITensorNetwork(to_callable(value), is; kwargs...)
-end
-
-function ITensorNetwork(
-        elt::Type, f, is::IndsNetwork; link_space = trivial_space(is), kwargs...
-    )
-    tn = ITensorNetwork(f, is; kwargs...)
-    for v in vertices(tn)
-        # TODO: Ideally we would use broadcasting, i.e. `elt.(tn[v])`,
-        # but that doesn't work right now on ITensors.
-        tn[v] = ITensors.convert_eltype(elt, tn[v])
-    end
-    return tn
-end
-
-function ITensorNetwork(
-        itensor_constructor::Function, is::IndsNetwork; link_space = trivial_space(is),
-        kwargs...
-    )
-    is = insert_linkinds(is; link_space)
-    tn = ITensorNetwork{vertextype(is)}(ITensor[])
-    for v in vertices(is)
-        add_vertex!(tn, v)
-    end
-    for e in edges(is)
-        add_edge!(tn, e)
-    end
-    for v in vertices(tn)
-        # TODO: Replace with `is[v]` once `getindex(::IndsNetwork, ...)` is smarter.
-        siteinds = get(is, v, Index[])
-        edges = [edgetype(is)(v, nv) for nv in neighbors(is, v)]
-        linkinds = map(e -> is[e], Indices(edges))
-        tensor_v = generic_state(itensor_constructor(v), (; siteinds, linkinds))
-        setindex_preserve_graph!(tn, tensor_v, v)
-    end
-    return tn
-end
diff --git a/src/solvers/eigsolve.jl b/src/solvers/eigsolve.jl
index dd741c61..2679147d 100644
--- a/src/solvers/eigsolve.jl
+++ b/src/solvers/eigsolve.jl
@@ -70,7 +70,7 @@ DMRG-like sweep algorithm on a `TreeTensorNetwork`.
 # Arguments
 
   - `operator`: The operator to diagonalize, typically a `TreeTensorNetwork` representing a
-    Hamiltonian constructed from an `OpSum` (e.g. via `ttn(opsum, sites)`).
+    Hamiltonian constructed from an `OpSum` (e.g. via `TreeTensorNetwork(opsum, sites)`).
   - `init_state`: Initial guess for the eigenvector as a `TreeTensorNetwork`.
   - `nsweeps`: Number of sweeps over the network.
   - `nsites=1`: Number of sites optimized simultaneously per local update step (1 or 2).
diff --git a/src/treetensornetworks/abstracttreetensornetwork.jl b/src/treetensornetworks/abstracttreetensornetwork.jl
index 8aa93cc4..85ebfa18 100644
--- a/src/treetensornetworks/abstracttreetensornetwork.jl
+++ b/src/treetensornetworks/abstracttreetensornetwork.jl
@@ -1,9 +1,11 @@
+using DataGraphs: underlying_graph
 using Graphs: has_vertex
-using ITensors: ITensors, @Algorithm_str, Algorithm, directsum, hasinds, permute, plev
+using ITensors:
+    ITensors, @Algorithm_str, Algorithm, ITensor, directsum, hasinds, permute, plev
 using IsApprox: IsApprox, Approx
 using NamedGraphs.GraphsExtensions: GraphsExtensions, a_star, edge_path, leaf_vertices,
-    post_order_dfs_edges, post_order_dfs_vertices
-using NamedGraphs: namedgraph_a_star, steiner_tree
+    post_order_dfs_edges, post_order_dfs_vertices, vertextype
+using NamedGraphs: NamedGraph, namedgraph_a_star, steiner_tree
 using TupleTools: TupleTools
 
 abstract type AbstractTreeTensorNetwork{V} <: AbstractITensorNetwork{V} end
@@ -271,8 +273,13 @@ function Base.:+(
         tns[j] = orthogonalize(tns[j], root_vertex)
     end
 
-    # Output state
-    tn = ttn(siteinds(tns[1]))
+    # Output state: empty TTN over the same graph as the inputs.
+    # Tensor data and link indices are filled in by the directsum loop below.
+    g_out = NamedGraph(underlying_graph(siteinds(tns[1])))
+    tensors_out = Dict{vertextype(g_out), ITensor}(
+        v => ITensor() for v in vertices(g_out)
+    )
+    tn = TreeTensorNetwork(ITensorNetwork(tensors_out, g_out))
 
     vs = post_order_dfs_vertices(tn, root_vertex)
     es = post_order_dfs_edges(tn, root_vertex)
diff --git a/src/treetensornetworks/opsum_to_ttn/opsum_to_ttn.jl b/src/treetensornetworks/opsum_to_ttn/opsum_to_ttn.jl
index 2f530858..0b607b3a 100644
--- a/src/treetensornetworks/opsum_to_ttn/opsum_to_ttn.jl
+++ b/src/treetensornetworks/opsum_to_ttn/opsum_to_ttn.jl
@@ -1,12 +1,13 @@
 #using FillArrays: OneElement
-#using DataGraphs: DataGraph
+using DataGraphs: underlying_graph
 using Graphs: degree, is_tree, rem_vertex!
 using ITensors.LazyApply: Prod, Sum, coefficient
 using ITensors.NDTensors: Block, blockdim, maxdim, nblocks, nnzblocks, truncate!
 using ITensors.Ops: Op, OpSum, argument, coefficient, name, params, site, terms, which_op
-using ITensors: flux, has_fermion_string, itensor, removeqns, space
+using ITensors: ITensor, flux, has_fermion_string, itensor, removeqns, space
 using NamedGraphs.GraphsExtensions:
-    GraphsExtensions, boundary_edges, degrees, is_leaf_vertex, vertex_path
+    GraphsExtensions, boundary_edges, degrees, is_leaf_vertex, vertex_path, vertextype
+using NamedGraphs: NamedGraph
 using StaticArrays: MVector
 
 #
@@ -289,7 +290,10 @@ function compress_ttn(
         )
     end
 
-    H = ttn(sites0)   # initialize TTN without the dummy indices added
+    # initialize TTN without the dummy indices added; tensors are filled in below
+    g0 = NamedGraph(underlying_graph(sites0))
+    tensors0 = Dict{vertextype(g0), ITensor}(v => ITensor() for v in vertices(g0))
+    H = TreeTensorNetwork(ITensorNetwork(tensors0, g0))
     function qnblock(i::Index, q::QN)
         for b in 2:(nblocks(i) - 1)
             flux(i, Block(b)) == q && return b
@@ -632,12 +636,12 @@ function sortmergeterms(os::OpSum{C}) where {C}
 end
 
 """
-    ttn(os::OpSum, sites::IndsNetwork{<:Index}; kwargs...)
-    ttn(eltype::Type{<:Number}, os::OpSum, sites::IndsNetwork{<:Index}; kwargs...)
+    TreeTensorNetwork(os::OpSum, sites::IndsNetwork{<:Index}; kwargs...)
+    TreeTensorNetwork(eltype::Type{<:Number}, os::OpSum, sites::IndsNetwork{<:Index}; kwargs...)
 
 Convert an OpSum object `os` to a TreeTensorNetwork, with indices given by `sites`.
 """
-function ttn(
+function TreeTensorNetwork(
         os::OpSum,
         sites::IndsNetwork;
         root_vertex = GraphsExtensions.default_root_vertex(sites),
@@ -654,16 +658,16 @@ function ttn(
 end
 
 function mpo(os::OpSum, external_inds::Vector; kwargs...)
-    return ttn(os, path_indsnetwork(external_inds); kwargs...)
+    return TreeTensorNetwork(os, path_indsnetwork(external_inds); kwargs...)
 end
 function mpo(os::OpSum, s::IndsNetwork; kwargs...)
     # TODO: Check it is a path graph.
-    return ttn(os, s; kwargs...)
+    return TreeTensorNetwork(os, s; kwargs...)
 end
 
 # Catch-all for leaf eltype specification
-function ttn(eltype::Type{<:Number}, os, sites::IndsNetwork; kwargs...)
-    return NDTensors.convert_scalartype(eltype, ttn(os, sites; kwargs...))
+function TreeTensorNetwork(eltype::Type{<:Number}, os, sites::IndsNetwork; kwargs...)
+    return NDTensors.convert_scalartype(eltype, TreeTensorNetwork(os, sites; kwargs...))
 end
 
 #
diff --git a/src/treetensornetworks/treetensornetwork.jl b/src/treetensornetworks/treetensornetwork.jl
index 48103384..7dc7cb25 100644
--- a/src/treetensornetworks/treetensornetwork.jl
+++ b/src/treetensornetworks/treetensornetwork.jl
@@ -1,9 +1,10 @@
+using DataGraphs: underlying_graph
 using Dictionaries: Indices
 using Graphs: path_graph
 using ITensors: ITensor
 using LinearAlgebra: factorize, normalize
 using NamedGraphs.GraphsExtensions: GraphsExtensions, vertextype
-using NamedGraphs: similar_graph
+using NamedGraphs: NamedGraph, similar_graph
 
 """
     TreeTensorNetwork{V} <: AbstractTreeTensorNetwork{V}
@@ -14,10 +15,10 @@ center of the network.
 
 `TTN` is an alias for `TreeTensorNetwork`.
 
-Use [`ttn`](@ref) or [`mps`](@ref) to construct instances, and [`orthogonalize`](@ref) to
-bring the network into a canonical gauge.
+Use the [`TreeTensorNetwork`](@ref) constructors to build instances, and
+[`orthogonalize`](@ref) to bring the network into a canonical gauge.
 
-See also: [`ITensorNetwork`](@ref), [`ttn`](@ref), [`mps`](@ref).
+See also: [`ITensorNetwork`](@ref).
 """
 struct TreeTensorNetwork{V} <: AbstractTreeTensorNetwork{V}
     tensornetwork::ITensorNetwork{V}
@@ -52,21 +53,25 @@ Throws an error if the underlying graph of `tn` is not a tree.
 ```jldoctest
 julia> using NamedGraphs.NamedGraphGenerators: named_comb_tree
 
+julia> using NamedGraphs: NamedGraph
+
 julia> using Graphs: vertices
 
+julia> using ITensors: ITensor
+
 julia> g = named_comb_tree((2, 2));
 
 julia> s = siteinds("S=1/2", g);
 
-julia> itn = ITensorNetwork(s; link_space = 2);
+julia> tensors = Dict(v => ITensor(s[v]...) for v in vertices(g));
 
-julia> root_vertex = first(vertices(itn));
+julia> itn = ITensorNetwork(tensors, NamedGraph(g));
 
-julia> ttn_state = TreeTensorNetwork(itn; ortho_region = [root_vertex]);
+julia> ttn_state = TreeTensorNetwork(itn; ortho_region = [first(vertices(itn))]);
 
 ```
 
-See also: [`ttn`](@ref), [`ITensorNetwork`](@ref), [`orthogonalize`](@ref).
+See also: [`ITensorNetwork`](@ref), [`orthogonalize`](@ref).
 """
 function TreeTensorNetwork(tn::ITensorNetwork; ortho_region = vertices(tn))
     return _TreeTensorNetwork(tn, ortho_region)
@@ -88,7 +93,7 @@ end
 Convert a `TreeTensorNetwork` to a plain `ITensorNetwork`, discarding orthogonality
 metadata. The returned network shares the same underlying tensor data.
 
-See also: [`TreeTensorNetwork`](@ref), [`ttn`](@ref).
+See also: [`TreeTensorNetwork`](@ref).
 """
 ITensorNetwork(tn::TTN) = copy(tn.tensornetwork)
 
@@ -119,79 +124,11 @@ end
 # set_ortho_region: low-level update of the ortho_region metadata only,
 # without any gauge transformations. To move the orthogonality center use orthogonalize.
 function set_ortho_region(tn::TTN, ortho_region)
-    return ttn(tn.tensornetwork; ortho_region)
+    return TreeTensorNetwork(tn.tensornetwork; ortho_region)
 end
 
 """
-    ttn(args...; ortho_region=nothing) -> TreeTensorNetwork
-
-Construct a `TreeTensorNetwork` (TTN) using the same interface as [`ITensorNetwork`](@ref).
-All positional and keyword arguments are forwarded to the `ITensorNetwork` constructor.
-
-If `ortho_region` is not specified, no particular gauge is assumed.
-Call [`orthogonalize`](@ref) to impose a gauge.
-
-# Example
-
-```jldoctest
-julia> using NamedGraphs.NamedGraphGenerators: named_comb_tree
-
-julia> g = named_comb_tree((2, 2));
-
-julia> s = siteinds("S=1/2", g);
-
-julia> psi = ttn(v -> "Up", s);
-
-```
-
-See also: [`mps`](@ref), [`TreeTensorNetwork`](@ref).
-"""
-function ttn(args...; ortho_region = nothing)
-    tn = ITensorNetwork(args...)
-    if isnothing(ortho_region)
-        ortho_region = vertices(tn)
-    end
-    return _TreeTensorNetwork(tn, ortho_region)
-end
-
-"""
-    mps(args...; ortho_region=nothing) -> TreeTensorNetwork
-
-Construct a matrix product state (MPS) as a `TreeTensorNetwork` on a 1D path graph.
-The interface is identical to [`ttn`](@ref) but is intended for 1D (chain) topologies.
-
-See also: [`ttn`](@ref).
-"""
-function mps(args...; ortho_region = nothing)
-    # TODO: Check it is a path graph.
-    tn = ITensorNetwork(args...)
-    if isnothing(ortho_region)
-        ortho_region = vertices(tn)
-    end
-    return _TreeTensorNetwork(tn, ortho_region)
-end
-
-"""
-    mps(f, is::Vector{<:Index}; kwargs...) -> TreeTensorNetwork
-
-Construct a matrix product state (MPS) from a function `f` and a flat vector of site
-indices `is`. The indices are arranged on a 1D path graph automatically.
-
-# Example
-
-```jldoctest
-julia> s = siteinds("S=1/2", 6);
-
-julia> psi = mps(v -> "Up", s);
-
-```
-"""
-function mps(f, is::Vector{<:Index}; kwargs...)
-    return mps(f, path_indsnetwork(is); kwargs...)
-end
-
-"""
-    ttn(a::ITensor, is::IndsNetwork; ortho_region=..., kwargs...) -> TreeTensorNetwork
+    TreeTensorNetwork(a::ITensor, is::IndsNetwork; ortho_region=..., kwargs...) -> TreeTensorNetwork
 
 Decompose a dense `ITensor` `a` into a `TreeTensorNetwork` with the tree structure
 described by the `IndsNetwork` `is`.
@@ -199,25 +136,8 @@ described by the `IndsNetwork` `is`.
 Successive QR/SVD factorizations are applied following a post-order DFS traversal from the
 root vertex, then the network is orthogonalized to `ortho_region` (defaults to the root).
 Extra `kwargs` (e.g. `cutoff`, `maxdim`) are forwarded to the factorization.
-
-# Example
-
-```jldoctest
-julia> using NamedGraphs.NamedGraphGenerators: named_comb_tree
-
-julia> using ITensors: ITensors
-
-julia> g = named_comb_tree((3, 1));
-
-julia> s = siteinds("S=1/2", g);
-
-julia> A = ITensors.random_itensor(only(s[(1, 1)]), only(s[(2, 1)]), only(s[(3, 1)]));
-
-julia> ttn_A = ttn(A, s);
-
-```
 """
-function ttn(
+function TreeTensorNetwork(
         a::ITensor,
         is::IndsNetwork;
         ortho_region = Indices([GraphsExtensions.default_root_vertex(is)]),
@@ -227,7 +147,16 @@ function ttn(
         @assert hasinds(a, is[v])
     end
     @assert ortho_region ⊆ vertices(is)
-    tn = ITensorNetwork(is)
+    is = insert_linkinds(is)
+    g = NamedGraph(underlying_graph(is))
+    ts = Dict{vertextype(g), ITensor}()
+    for v in vertices(g)
+        site_inds = get(is, v, Index[])
+        edges_v = [edgetype(is)(v, nv) for nv in neighbors(is, v)]
+        link_inds = reduce(vcat, (is[e] for e in edges_v); init = Index[])
+        ts[v] = ITensor(site_inds..., link_inds...)
+    end
+    tn = ITensorNetwork(ts, g)
     ortho_center = first(ortho_region)
     for e in post_order_dfs_edges(tn, ortho_center)
         left_inds = setdiff(inds(a), inds(tn[dst(e)]))
@@ -237,6 +166,6 @@ function ttn(
         a = a_r
     end
     tn[ortho_center] = a
-    ttn_a = ttn(tn)
+    ttn_a = TreeTensorNetwork(tn)
     return orthogonalize(ttn_a, ortho_center)
 end
diff --git a/test/Project.toml b/test/Project.toml
index 708ec90f..11505d90 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -5,7 +5,6 @@ Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"
 DataGraphs = "b5a273c3-7e6c-41f6-98bd-8d7f1525a36a"
 Dictionaries = "85a47980-9c8c-11e8-2b9f-f7ca1fa99fb4"
-Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 Glob = "c27321d9-0574-5035-807b-f59d2c89b15c"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
 GraphsFlows = "06909019-6f44-4949-96fc-b9d9aaa02889"
@@ -40,12 +39,11 @@ Aqua = "0.8.11"
 Compat = "4.16"
 DataGraphs = "0.4"
 Dictionaries = "0.4.4"
-Distributions = "0.25.118"
 Glob = "1.3.1"
 Graphs = "1.12"
 GraphsFlows = "0.1.1"
 ITensorMPS = "0.3.6, 0.4"
-ITensorNetworks = "0.21"
+ITensorNetworks = "0.22"
 ITensorPkgSkeleton = "0.3.42"
 ITensors = "0.7, 0.8, 0.9"
 KrylovKit = "0.8, 0.9, 0.10"
diff --git a/test/solvers/test_applyexp.jl b/test/solvers/test_applyexp.jl
index 81f0e0e4..faa4e1db 100644
--- a/test/solvers/test_applyexp.jl
+++ b/test/solvers/test_applyexp.jl
@@ -1,12 +1,13 @@
 using Graphs: add_edge!, add_vertex!, vertices
 using ITensorNetworks:
-    ITensorNetworks, applyexp, dmrg, maxlinkdim, siteinds, time_evolve, ttn
+    ITensorNetworks, TreeTensorNetwork, applyexp, dmrg, maxlinkdim, siteinds, time_evolve
 using ITensors
 using ITensors.Ops: OpSum
 using NamedGraphs.NamedGraphGenerators: named_path_graph
 using NamedGraphs: NamedGraph
 using TensorOperations: TensorOperations
 using Test: @test, @testset #For contraction order finding
+include(joinpath(@__DIR__, "..", "utils.jl"))
 
 function chain_plus_ancilla(; nchain)
     g = NamedGraph()
@@ -38,14 +39,14 @@ end
             h += 1 / 2, "S+", j, "S-", j + 1
             h += 1 / 2, "S-", j, "S+", j + 1
         end
-        H = ttn(h, sites)
+        H = TreeTensorNetwork(h, sites)
 
         # Make initial product state
         state = Dict{Int, String}()
         for (j, v) in enumerate(vertices(sites))
             state[v] = iseven(j) ? "Up" : "Dn"
         end
-        psi0 = ttn(state, sites)
+        psi0 = TreeTensorNetwork(productstate(state, sites))
 
         cutoff = 1.0e-10
         maxdim = 100
@@ -89,14 +90,14 @@ end
             h += 1 / 2, "S+", j, "S-", j + 1
             h += 1 / 2, "S-", j, "S+", j + 1
         end
-        H = ttn(h, sites)
+        H = TreeTensorNetwork(h, sites)
 
         # Initial product state
         state = Dict{Int, String}()
         for (j, v) in enumerate(vertices(sites))
             state[v] = iseven(j) ? "Up" : "Dn"
         end
-        psi0 = ttn(state, sites)
+        psi0 = TreeTensorNetwork(productstate(state, sites))
 
         nsites = 2
         factorize_kwargs = (; cutoff = 1.0e-8, maxdim = 100)
diff --git a/test/solvers/test_eigsolve.jl b/test/solvers/test_eigsolve.jl
index 6be2e4bd..5805fdf1 100644
--- a/test/solvers/test_eigsolve.jl
+++ b/test/solvers/test_eigsolve.jl
@@ -1,11 +1,12 @@
 using Graphs: dst, edges, src, vertices
-using ITensorNetworks: SweepIterator, dmrg, siteinds, ttn
+using ITensorNetworks: SweepIterator, TreeTensorNetwork, dmrg, siteinds
 using ITensors
 using ITensors.Ops: OpSum
 using Suppressor: @capture_out
 using TensorOperations: TensorOperations
 using Test: @test, @testset
 
+include(joinpath(@__DIR__, "..", "utils.jl"))
 include("utilities/simple_ed_methods.jl")
 include("utilities/tree_graphs.jl")
 
@@ -24,14 +25,14 @@ include("utilities/tree_graphs.jl")
         h += 1 / 2, "S+", i, "S-", j
         h += 1 / 2, "S-", i, "S+", j
     end
-    H = ttn(h, sites)
+    H = TreeTensorNetwork(h, sites)
 
     # Make initial product state
     state = Dict{Tuple{Int, Int}, String}()
     for (j, v) in enumerate(vertices(sites))
         state[v] = iseven(j) ? "Up" : "Dn"
     end
-    psi0 = ttn(state, sites)
+    psi0 = TreeTensorNetwork(productstate(state, sites))
 
     (outputlevel >= 1) && println("Computing exact ground state")
     Ex, psix = ed_ground_state(H, psi0)
diff --git a/test/test_additensornetworks.jl b/test/test_additensornetworks.jl
index b9676862..7503fd1f 100644
--- a/test/test_additensornetworks.jl
+++ b/test/test_additensornetworks.jl
@@ -1,5 +1,6 @@
 using Graphs: rem_edge!, vertices
-using ITensorNetworks: ITensorNetwork, inner_network, orthogonalize, siteinds, truncate, ttn
+using ITensorNetworks:
+    ITensorNetwork, TreeTensorNetwork, inner_network, orthogonalize, siteinds, truncate
 using ITensors: ITensors, apply, inner, op, scalar
 using LinearAlgebra: norm_sqr
 using NamedGraphs.NamedGraphGenerators: named_comb_tree, named_grid
@@ -11,8 +12,8 @@ include("utils.jl")
 @testset "add_itensornetworks" begin
     g = named_grid((2, 2))
     s = siteinds("S=1/2", g)
-    ψ1 = ITensorNetwork(v -> "↑", s)
-    ψ2 = ITensorNetwork(v -> "↓", s)
+    ψ1 = productstate(v -> "↑", s)
+    ψ2 = productstate(v -> "↓", s)
 
     ψ_GHZ = ψ1 + ψ2
 
@@ -71,13 +72,13 @@ end
     for (j, v) in enumerate(verts)
         state1[v] = isodd(j) ? "Up" : "Dn"
     end
-    ψ1 = ttn(state1, sites)
+    ψ1 = TreeTensorNetwork(productstate(state1, sites))
 
     state2 = Dict{Tuple, String}()
     for (j, v) in enumerate(vertices(g))
         state2[v] = isodd(j) ? "Dn" : "Up"
     end
-    ψ2 = ttn(state2, sites)
+    ψ2 = TreeTensorNetwork(productstate(state2, sites))
 
     ϕ = ψ1 + ψ2
 
diff --git a/test/test_expect.jl b/test/test_expect.jl
index 73a0891b..98040d98 100644
--- a/test/test_expect.jl
+++ b/test/test_expect.jl
@@ -43,7 +43,7 @@ include("utils.jl")
     g = named_grid((L, L))
     s = siteinds("S=1/2", g; conserve_qns = true)
 
-    ψ = ITensorNetwork(v -> isodd(sum(v)) ? "↑" : "↓", s)
+    ψ = productstate(v -> isodd(sum(v)) ? "↑" : "↓", s)
 
     sz_bp = expect(ψ, "Sz"; alg = "bp", cache_update_kwargs = (; maxiter = 20))
     sz_exact = expect(ψ, "Sz"; alg = "exact")
diff --git a/test/test_inner.jl b/test/test_inner.jl
index 0e0a796c..81df35a5 100644
--- a/test/test_inner.jl
+++ b/test/test_inner.jl
@@ -1,6 +1,6 @@
 using Graphs: SimpleGraph, uniform_tree
-using ITensorNetworks: ITensorNetwork, inner, inner_network, loginner, logscalar, scalar,
-    siteinds, ttn, underlying_graph
+using ITensorNetworks: ITensorNetwork, TreeTensorNetwork, inner, inner_network, loginner,
+    logscalar, scalar, siteinds, underlying_graph
 using ITensors: dag
 using NamedGraphs: NamedGraph
 using SplitApplyCombine: group
@@ -38,7 +38,7 @@ using .ModelHamiltonians: heisenberg
     @test xy_scalar ≈ xy_scalar_logbp
 
     #test contraction of three layers for expectation values
-    A = ITensorNetwork(ttn(heisenberg(g), s))
+    A = ITensorNetwork(TreeTensorNetwork(heisenberg(g), s))
     xAy_scalar = inner(x, A, y; alg = "exact")
     xAy_scalar_bp = inner(x, A, y; alg = "bp")
     xAy_scalar_logbp = exp(loginner(x, A, y; alg = "bp"))
diff --git a/test/test_itensornetwork.jl b/test/test_itensornetwork.jl
index 7f1b6926..05c6dda1 100644
--- a/test/test_itensornetwork.jl
+++ b/test/test_itensornetwork.jl
@@ -1,11 +1,10 @@
 using Dictionaries: Dictionary
-using Distributions: Uniform
 using Graphs: degree, dijkstra_shortest_paths, edges, grid, has_vertex, ne, neighbors, nv,
     rem_vertex!, vertices, weights
 using GraphsFlows: GraphsFlows
 using ITensorNetworks: ITensorNetworks, ITensorNetwork, IndsNetwork, contraction_sequence,
     inner_network, linkinds, norm_sqr, norm_sqr_network, orthogonalize, siteinds,
-    tree_orthogonalize, ttn, ⊗
+    tree_orthogonalize, ⊗
 using ITensors.NDTensors: NDTensors, dim
 using ITensors: ITensors, ITensor, Index, Op, commonind, commoninds, contract, dag, hasinds,
     inds, inner, itensor, onehot, order, prime, random_itensor, scalartype, sim
@@ -28,7 +27,8 @@ const elts = (Float32, Float64, Complex{Float32}, Complex{Float64})
         @test nv(s) == 4
         @test ne(s) == 3
         @test neighbors(s, (2,)) == [(1,), (3,)]
-        tn = ITensorNetwork(s; link_space = 2)
+        rng = StableRNG(1234)
+        tn = random_tensornetwork(rng, s; link_space = 2)
         @test nv(tn) == 4
         @test ne(tn) == 3
         @test tn isa ITensorNetwork
@@ -54,10 +54,6 @@ const elts = (Float32, Float64, Complex{Float32}, Complex{Float64})
         @test_broken tn[1:2] isa ITensorNetwork
         # TODO: Support this syntax, maybe rename `subgraph`.
         @test_broken induced_subgraph(tn, [(1,), (2,)]) isa ITensorNetwork
-        rng = StableRNG(1234)
-        for v in vertices(tn)
-            tn[v] = randn!(rng, tn[v])
-        end
         tn′ = sim(dag(tn); sites = [])
 
         @test tn′ isa ITensorNetwork
@@ -69,7 +65,7 @@ const elts = (Float32, Float64, Complex{Float32}, Complex{Float64})
         @test inner_res isa Float64
 
         # test that by default vertices are linked by bond-dimension 1 index
-        tn = ITensorNetwork(s)
+        tn = random_tensornetwork(rng, s)
         @test isone(ITensors.dim(commonind(tn[(1,)], tn[(2,)])))
     end
 
@@ -92,7 +88,7 @@ const elts = (Float32, Float64, Complex{Float32}, Complex{Float64})
         dims = (2, 2)
         g = named_grid(dims)
         s = siteinds("S=1/2", g)
-        ψ = ITensorNetwork(v -> "↑", s)
+        ψ = productstate(v -> "↑", s)
         tn = disjoint_union("bra" => ψ, "ket" => prime(dag(ψ); sites = []))
         tn_2 = contract(tn, ((1, 2), "ket") => ((1, 2), "bra"))
         @test !has_vertex(tn_2, ((1, 2), "ket"))
@@ -103,7 +99,7 @@ const elts = (Float32, Float64, Complex{Float32}, Complex{Float64})
         dims = (2, 2)
         g = named_grid(dims)
         s = siteinds("S=1/2", g)
-        ψ = ITensorNetwork(v -> "↑", s)
+        ψ = productstate(v -> "↑", s)
         rem_vertex!(ψ, (1, 2))
         tn = norm_sqr_network(ψ)
         @test !has_vertex(tn, ((1, 2), "bra"))
@@ -126,9 +122,7 @@ const elts = (Float32, Float64, Complex{Float32}, Complex{Float64})
             )
 
         rng = StableRNG(1234)
-        ψ = ITensorNetwork(g; kwargs...) do v
-            return inds -> itensor(randn(rng, elt, dim.(inds)...), inds)
-        end
+        ψ = random_tensornetwork(rng, elt, g; kwargs...)
         @test eltype(ψ[first(vertices(ψ))]) == elt
 
         ψc = conj(ψ)
@@ -141,97 +135,11 @@ const elts = (Float32, Float64, Complex{Float32}, Complex{Float64})
             @test ψd[v] == dag(ψ[v])
         end
 
-        rng = StableRNG(1234)
-        ψ = ITensorNetwork(g; kwargs...) do v
-            return inds -> itensor(randn(rng, dim.(inds)...), inds)
-        end
-        @test eltype(ψ[first(vertices(ψ))]) == Float64
-        rng = StableRNG(1234)
-        ψ = random_tensornetwork(rng, elt, g; kwargs...)
-        @test eltype(ψ[first(vertices(ψ))]) == elt
         rng = StableRNG(1234)
         ψ = random_tensornetwork(rng, g; kwargs...)
         @test eltype(ψ[first(vertices(ψ))]) == Float64
-        ψ = ITensorNetwork(elt, undef, g; kwargs...)
-        @test eltype(ψ[first(vertices(ψ))]) == elt
-        ψ = ITensorNetwork(undef, g)
-        @test eltype(ψ[first(vertices(ψ))]) == Float64
     end
 
-    @testset "Product state constructors" for elt in elts
-        dims = (2, 2)
-        g = named_comb_tree(dims)
-        s = siteinds("S=1/2", g)
-        state1 = ["↑" "↓"; "↓" "↑"]
-        state2 = reshape([[1, 0], [0, 1], [0, 1], [1, 0]], 2, 2)
-        each_args = (;
-            ferro = (
-                ("↑",),
-                (elt, "↑"),
-                (Returns(i -> ITensor([1, 0], i)),),
-                (elt, Returns(i -> ITensor([1, 0], i))),
-                (Returns([1, 0]),),
-                (elt, Returns([1, 0])),
-            ),
-            antiferro = (
-                (state1,),
-                (elt, state1),
-                (Dict(CartesianIndices(dims) .=> state1),),
-                (elt, Dict(CartesianIndices(dims) .=> state1)),
-                (Dict(Tuple.(CartesianIndices(dims)) .=> state1),),
-                (elt, Dict(Tuple.(CartesianIndices(dims)) .=> state1)),
-                (Dictionary(CartesianIndices(dims), state1),),
-                (elt, Dictionary(CartesianIndices(dims), state1)),
-                (Dictionary(Tuple.(CartesianIndices(dims)), state1),),
-                (elt, Dictionary(Tuple.(CartesianIndices(dims)), state1)),
-                (state2,),
-                (elt, state2),
-                (Dict(CartesianIndices(dims) .=> state2),),
-                (elt, Dict(CartesianIndices(dims) .=> state2)),
-                (Dict(Tuple.(CartesianIndices(dims)) .=> state2),),
-                (elt, Dict(Tuple.(CartesianIndices(dims)) .=> state2)),
-                (Dictionary(CartesianIndices(dims), state2),),
-                (elt, Dictionary(CartesianIndices(dims), state2)),
-                (Dictionary(Tuple.(CartesianIndices(dims)), state2),),
-                (elt, Dictionary(Tuple.(CartesianIndices(dims)), state2)),
-            ),
-        )
-        for pattern in keys(each_args)
-            for args in each_args[pattern]
-                x = ITensorNetwork(args..., s)
-                if first(args) === elt
-                    @test scalartype(x) === elt
-                else
-                    @test scalartype(x) === Float64
-                end
-                for v in vertices(x)
-                    xᵛ = x[v]
-                    @test degree(x, v) + 1 == ndims(xᵛ)
-                    sᵛ = only(siteinds(x, v))
-                    for w in neighbors(x, v)
-                        lʷ = only(linkinds(x, v => w))
-                        @test dim(lʷ) == 1
-                        xᵛ *= onehot(lʷ => 1)
-                    end
-                    @test ndims(xᵛ) == 1
-                    a = if pattern == :ferro
-                        [1, 0]
-                    elseif pattern == :antiferro
-                        iseven(sum(v)) ? [1, 0] : [0, 1]
-                    end
-                    @test xᵛ == ITensor(a, sᵛ)
-                end
-            end
-        end
-    end
-    @testset "random_tensornetwork with custom distributions" begin
-        distribution = Uniform(-1.0, 1.0)
-        rng = StableRNG(1234)
-        tn = random_tensornetwork(rng, distribution, named_grid(4); link_space = 2)
-        # Note: distributions in package `Distributions` currently doesn't support customized
-        # eltype, and all elements have type `Float64`
-        @test eltype(tn[first(vertices(tn))]) == Float64
-    end
     @testset "orthogonalize" begin
         rng = StableRNG(1234)
         tn = random_tensornetwork(rng, named_grid(4); link_space = 2)
@@ -256,7 +164,8 @@ const elts = (Float32, Float64, Complex{Float32}, Complex{Float64})
     end
 
     @testset "dijkstra_shortest_paths" begin
-        tn = ITensorNetwork(named_grid(4); link_space = 2)
+        rng = StableRNG(1234)
+        tn = random_tensornetwork(rng, named_grid(4); link_space = 2)
         paths = dijkstra_shortest_paths(tn, [1])
         @test paths.dists == Dictionary([0, 1, 2, 3])
         @test paths.parents == Dictionary([1, 1, 2, 3])
@@ -264,7 +173,8 @@ const elts = (Float32, Float64, Complex{Float32}, Complex{Float64})
     end
 
     @testset "mincut" begin
-        tn = ITensorNetwork(named_grid(4); link_space = 3)
+        rng = StableRNG(1234)
+        tn = random_tensornetwork(rng, named_grid(4); link_space = 3)
         w = weights(tn)
         @test w isa Dictionary{Tuple{Int, Int}, Float64}
         @test length(w) ≈ ne(tn)
@@ -286,7 +196,8 @@ const elts = (Float32, Float64, Complex{Float32}, Complex{Float64})
         dims = (2, 2)
         g = named_grid(dims)
         s = siteinds("S=1/2", g)
-        ψ = ITensorNetwork(s; link_space = 2)
+        rng = StableRNG(1234)
+        ψ = random_tensornetwork(rng, s; link_space = 2)
 
         e = (1, 1) => (2, 1)
         uie = setdiff(inds(ψ[1, 1]), inds(ψ[2, 1]))
@@ -321,14 +232,14 @@ const elts = (Float32, Float64, Complex{Float32}, Complex{Float64})
         g = named_grid(dims)
         s = siteinds("S=1/2", g)
         state_map(v::Tuple) = iseven(sum(isodd.(v))) ? "↑" : "↓"
-        ψ = ITensorNetwork(state_map, s)
+        ψ = productstate(state_map, s)
         t = ψ[2, 2]
         si = only(siteinds(ψ, (2, 2)))
         bi = map(e -> only(linkinds(ψ, e)), incident_edges(ψ, (2, 2)))
         @test eltype(t) == Float64
         @test abs(t[si => "↑", [b => end for b in bi]...]) == 1.0 # insert_links introduces extra signs through factorization...
         @test t[si => "↓", [b => end for b in bi]...] == 0.0
-        ϕ = ITensorNetwork(elt, state_map, s)
+        ϕ = productstate(elt, state_map, s)
         t = ϕ[2, 2]
         si = only(siteinds(ϕ, (2, 2)))
         bi = map(e -> only(linkinds(ϕ, e)), incident_edges(ϕ, (2, 2)))
diff --git a/test/test_opsum_to_ttn.jl b/test/test_opsum_to_ttn.jl
index f54076e8..ae8eb67a 100644
--- a/test/test_opsum_to_ttn.jl
+++ b/test/test_opsum_to_ttn.jl
@@ -1,5 +1,5 @@
 using Graphs: vertices
-using ITensorNetworks: ITensorNetworks, OpSum, siteinds, ttn
+using ITensorNetworks: ITensorNetworks, OpSum, TreeTensorNetwork, siteinds
 using ITensors.NDTensors: with_auto_fermion
 using NamedGraphs.NamedGraphGenerators: named_grid
 using Test: @test, @testset
@@ -15,9 +15,9 @@ using Test: @test, @testset
             os1 += 1.0, "Sx", (1, 1)
             os2 = OpSum()
             os2 += 1.0, "Sy", (1, 1)
-            H1 = ttn(os1, s)
-            H2 = ttn(os2, s)
-            H3 = ttn(os1 + os2, s)
+            H1 = TreeTensorNetwork(os1, s)
+            H2 = TreeTensorNetwork(os2, s)
+            H3 = TreeTensorNetwork(os1 + os2, s)
 
             @test H1 + H2 ≈ H3 rtol = 1.0e-6
         end
diff --git a/test/test_opsum_to_ttn_mpo_cross_check.jl b/test/test_opsum_to_ttn_mpo_cross_check.jl
index a7b16143..a0509f44 100644
--- a/test/test_opsum_to_ttn_mpo_cross_check.jl
+++ b/test/test_opsum_to_ttn_mpo_cross_check.jl
@@ -3,7 +3,7 @@ using Dictionaries: Dictionary
 using Graphs: add_edge!, add_vertex!, rem_edge!, vertices
 using ITensorMPS: ITensorMPS
 using ITensorNetworks.ITensorsExtensions: replace_vertices
-using ITensorNetworks: ITensorNetworks, siteinds, ttn
+using ITensorNetworks: ITensorNetworks, TreeTensorNetwork, siteinds
 using ITensors.NDTensors: matrix, with_auto_fermion
 using ITensors: @disable_warn_order, ITensor, Index, combinedind, combiner, contract, dag,
     inds, removeqns
@@ -53,7 +53,7 @@ end
 
         @testset "Svd approach" for root_vertex in leaf_vertices(is)
             # get TTN Hamiltonian directly
-            Hsvd = ttn(H, is; root_vertex, cutoff = 1.0e-10)
+            Hsvd = TreeTensorNetwork(H, is; root_vertex, cutoff = 1.0e-10)
             # get corresponding MPO Hamiltonian
             Hline = ITensorMPS.MPO(replace_vertices(v -> vmap[v], H), sites)
             # compare resulting dense Hamiltonians
@@ -63,7 +63,7 @@ end
             end
             @test Tttno ≈ Tmpo rtol = 1.0e-6
 
-            Hsvd_lr = ttn(Hlr, is; root_vertex, cutoff = 1.0e-10)
+            Hsvd_lr = TreeTensorNetwork(Hlr, is; root_vertex, cutoff = 1.0e-10)
             Hline_lr = ITensorMPS.MPO(replace_vertices(v -> vmap[v], Hlr), sites)
             @disable_warn_order begin
                 Tttno_lr = prod(Hline_lr)
@@ -102,7 +102,7 @@ end
 
         @testset "Svd approach" for root_vertex in leaf_vertices(is)
             # get TTN Hamiltonian directly
-            Hsvd = ttn(H, is; root_vertex, cutoff = 1.0e-10)
+            Hsvd = TreeTensorNetwork(H, is; root_vertex, cutoff = 1.0e-10)
             # get corresponding MPO Hamiltonian
             Hline = ITensorMPS.MPO(replace_vertices(v -> vmap[v], H), sites)
             # compare resulting sparse Hamiltonians
@@ -113,7 +113,7 @@ end
             end
             @test Tttno ≈ Tmpo rtol = 1.0e-6
 
-            Hsvd_lr = ttn(Hlr, is; root_vertex, cutoff = 1.0e-10)
+            Hsvd_lr = TreeTensorNetwork(Hlr, is; root_vertex, cutoff = 1.0e-10)
             Hline_lr = ITensorMPS.MPO(replace_vertices(v -> vmap[v], Hlr), sites)
             @disable_warn_order begin
                 Tttno_lr = prod(Hline_lr)
@@ -142,7 +142,7 @@ end
 
             @testset "Svd approach" for root_vertex in leaf_vertices(is)
                 # get TTN Hamiltonian directly
-                Hsvd = ttn(H, is; root_vertex, cutoff = 1.0e-10)
+                Hsvd = TreeTensorNetwork(H, is; root_vertex, cutoff = 1.0e-10)
                 # get corresponding MPO Hamiltonian
                 sites =
                     [only(is[v]) for v in reverse(post_order_dfs_vertices(c, root_vertex))]
@@ -209,7 +209,7 @@ end
 
         @testset "Svd approach" for root_vertex in leaf_vertices(is)
             # get TTN Hamiltonian directly
-            Hsvd = ttn(H, is_missing_site; root_vertex, cutoff = 1.0e-10)
+            Hsvd = TreeTensorNetwork(H, is_missing_site; root_vertex, cutoff = 1.0e-10)
             # get corresponding MPO Hamiltonian
             Hline = ITensorMPS.MPO(replace_vertices(v -> vmap[v], H), sites)
 
@@ -220,7 +220,7 @@ end
             end
             @test Tttno ≈ Tmpo rtol = 1.0e-6
 
-            Hsvd_lr = ttn(Hlr, is_missing_site; root_vertex, cutoff = 1.0e-10)
+            Hsvd_lr = TreeTensorNetwork(Hlr, is_missing_site; root_vertex, cutoff = 1.0e-10)
             Hline_lr = ITensorMPS.MPO(replace_vertices(v -> vmap[v], Hlr), sites)
             @disable_warn_order begin
                 Tttno_lr = prod(Hline_lr)
diff --git a/test/test_tebd.jl b/test/test_tebd.jl
index 8ff87991..17993cb2 100644
--- a/test/test_tebd.jl
+++ b/test/test_tebd.jl
@@ -27,7 +27,7 @@ ITensors.disable_warn_order()
     β = 2.0
     Δβ = 0.2
 
-    ψ_init = ITensorNetwork(v -> "↑", s)
+    ψ_init = productstate(v -> "↑", s)
     #E0 = expect(ℋ, ψ_init)
     ψ = tebd(
         group_terms(ℋ, g),
diff --git a/test/test_ttn_expect.jl b/test/test_ttn_expect.jl
index e9fabd7d..2ba474cb 100644
--- a/test/test_ttn_expect.jl
+++ b/test/test_ttn_expect.jl
@@ -1,9 +1,10 @@
 using Graphs: vertices
-using ITensorNetworks: expect, siteinds, ttn
+using ITensorNetworks: TreeTensorNetwork, expect, siteinds
 using LinearAlgebra: norm
 using NamedGraphs.NamedGraphGenerators: named_comb_tree
 using StableRNGs: StableRNG
 using Test: @test, @testset
+include("utils.jl")
 @testset "TTN expect" begin
     tooth_lengths = fill(2, 2)
     c = named_comb_tree(tooth_lengths)
@@ -15,7 +16,7 @@ using Test: @test, @testset
         magnetization[v] = isodd(i) ? 0.5 : -0.5
     end
     states = v -> d[v]
-    state = ttn(states, s)
+    state = TreeTensorNetwork(productstate(states, s))
     res = expect("Sz", state)
     @test all([isapprox(res[v], magnetization[v]; atol = 1.0e-8) for v in vertices(s)])
 end
diff --git a/test/test_ttn_position.jl b/test/test_ttn_position.jl
index 5be35608..f2e0e3ab 100644
--- a/test/test_ttn_position.jl
+++ b/test/test_ttn_position.jl
@@ -1,10 +1,11 @@
 using Dictionaries: Dictionary, Indices
 using Graphs: vertices
-using ITensorNetworks: ITensorNetwork, ProjTTN, environments, position, siteinds, ttn
+using ITensorNetworks:
+    ITensorNetwork, ProjTTN, TreeTensorNetwork, environments, position, siteinds
 using ITensors.NDTensors: with_auto_fermion
 using ITensors: ITensor
 using NamedGraphs.NamedGraphGenerators: named_comb_tree, named_path_graph
-using NamedGraphs: NamedEdge
+using NamedGraphs: NamedEdge, NamedGraph
 using Test: @test, @testset
 include("utils.jl")
 using .ModelHamiltonians: ModelHamiltonians
@@ -24,14 +25,14 @@ using .ModelHamiltonians: ModelHamiltonians
 
         os = ModelHamiltonians.heisenberg(c)
 
-        H = ttn(os, s)
+        H = TreeTensorNetwork(os, s)
 
         d = Dict()
         for (i, v) in enumerate(vertices(s))
             d[v] = isodd(i) ? "Up" : "Dn"
         end
         states = v -> d[v]
-        psi = ttn(states, s)
+        psi = TreeTensorNetwork(productstate(states, s))
 
         # actual test, verifies that position is out of place
         vs = collect(vertices(s))
@@ -46,8 +47,9 @@ using .ModelHamiltonians: ModelHamiltonians
 end
 @testset "ProjTTN construction regression test" begin
     pos = Indices{Tuple{String, Int}}()
-    g = named_path_graph(2)
-    operator = ttn(ITensorNetwork{Any}(g))
+    g = NamedGraph{Any}(named_path_graph(2))
+    tensors = Dict{Any, ITensor}(v => ITensor() for v in vertices(g))
+    operator = TreeTensorNetwork(ITensorNetwork{Any}(tensors, g))
     environments = Dictionary{NamedEdge{Any}, ITensor}()
     @test ProjTTN(pos, operator, environments) isa ProjTTN{Any, Indices{Any}}
 end
diff --git a/test/test_ttno.jl b/test/test_ttno.jl
index b4e8b55f..c876a3d9 100644
--- a/test/test_ttno.jl
+++ b/test/test_ttno.jl
@@ -1,5 +1,5 @@
 using Graphs: vertices
-using ITensorNetworks: contract, ortho_region, siteinds, ttn, union_all_inds
+using ITensorNetworks: TreeTensorNetwork, contract, ortho_region, siteinds, union_all_inds
 using ITensors: @disable_warn_order, prime, random_itensor
 using LinearAlgebra: norm
 using NamedGraphs.NamedGraphGenerators: named_comb_tree
@@ -25,7 +25,7 @@ using Test: @test, @testset
         rng = StableRNG(1234)
         O = random_itensor(rng, sites_o...)
         # dense TTN constructor from IndsNetwork
-        @disable_warn_order o1 = ttn(O, is_isp; cutoff)
+        @disable_warn_order o1 = TreeTensorNetwork(O, is_isp; cutoff)
         root_vertex = only(ortho_region(o1))
         @disable_warn_order begin
             O1 = contract(o1, root_vertex)
diff --git a/test/test_ttns.jl b/test/test_ttns.jl
index 1c045db9..64fac579 100644
--- a/test/test_ttns.jl
+++ b/test/test_ttns.jl
@@ -1,13 +1,14 @@
 using DataGraphs: vertex_data
 using Graphs: vertices
 using ITensorNetworks:
-    ITensorNetwork, TreeTensorNetwork, contract, ortho_region, orthogonalize, siteinds, ttn
+    ITensorNetwork, TreeTensorNetwork, contract, ortho_region, orthogonalize, siteinds
 using ITensors: @disable_warn_order, random_itensor
 using LinearAlgebra: norm
 using NamedGraphs.NamedGraphGenerators: named_comb_tree
 using Random: shuffle
 using StableRNGs: StableRNG
 using Test: @test, @testset
+include("utils.jl")
 @testset "TTN Basics" begin
     # random comb tree
     rng = StableRNG(1234)
@@ -25,7 +26,7 @@ using Test: @test, @testset
         rng = StableRNG(1234)
         S = random_itensor(rng, vertex_data(is)...)
         # dense TTN constructor from IndsNetwork
-        @disable_warn_order s1 = ttn(S, is; cutoff)
+        @disable_warn_order s1 = TreeTensorNetwork(S, is; cutoff)
         root_vertex = only(ortho_region(s1))
         @disable_warn_order begin
             S1 = contract(s1, root_vertex)
@@ -37,8 +38,9 @@ using Test: @test, @testset
         g = named_comb_tree((3, 2))
         sites = siteinds("S=1/2", g)
 
-        psi = ttn(sites)  # zero-initialised
-        psi = ttn(v -> "Up", sites)  # product state
+        rng = StableRNG(1234)
+        psi = TreeTensorNetwork(random_tensornetwork(rng, sites))  # random
+        psi = TreeTensorNetwork(productstate(v -> "Up", sites))  # product state
 
         itn = ITensorNetwork(psi)  # TTN → ITensorNetwork
         @test vertex_data(itn) == vertex_data(psi.tensornetwork)
@@ -50,8 +52,9 @@ using Test: @test, @testset
         g = named_comb_tree((3, 2))
         sites = siteinds("S=1/2", g)
 
-        psi = ttn(sites)  # zero-initialised
-        psi = ttn(v -> "Up", sites)  # product state
+        rng = StableRNG(1234)
+        psi = TreeTensorNetwork(random_tensornetwork(rng, sites))  # random
+        psi = TreeTensorNetwork(productstate(v -> "Up", sites))  # product state
 
         v1 = collect(vertices(psi))[1]
         v2 = collect(vertices(psi))[2]
diff --git a/test/utils.jl b/test/utils.jl
index 69a54020..ffcec865 100644
--- a/test/utils.jl
+++ b/test/utils.jl
@@ -3,120 +3,116 @@
 # `test_*.jl`). Each test file that needs these does `include("utils.jl")`
 # inside its gensym module.
 
-using DataGraphs: IsUnderlyingGraph, vertex_data
-using Distributions: Distribution
-using Graphs: edges, vertices
-using ITensorNetworks: ITensorNetworks, ITensorNetwork, IndsNetwork
+using DataGraphs: underlying_graph, vertex_data
+using Graphs: AbstractGraph, dst, edges, src, vertices
+using ITensorNetworks.ITensorsExtensions: trivial_space
+using ITensorNetworks: ITensorNetwork, IndsNetwork
 using ITensors.NDTensors: dim
-using ITensors: ITensor, Index, itensor
-using LinearAlgebra: normalize
+using ITensors: ITensors, ITensor, Index, dag, itensor, onehot
 using NamedGraphs.GraphsExtensions: incident_edges
+using NamedGraphs: NamedGraph
 using Random: Random, AbstractRNG
-using SimpleTraits: SimpleTraits, @traitfn
-
-# Build an `ITensorNetwork` on the graph backing the inds network `s` from a
-# user-supplied per-vertex tensor builder. Each vertex `v` is given a tensor
-# whose indices are `s[v]` (the site indices) followed by one fresh `Index` of
-# dimension `link_space` for each edge of the graph incident to `v`. The shared
-# link index is used on both endpoints so the resulting network's edges agree
-# with `edges(s)`.
-function _random_tensornetwork(s::IndsNetwork, build_tensor; link_space = 1)
-    l = Dict(e => Index(link_space, "Link") for e in edges(s))
-    l = merge(l, Dict(reverse(e) => l[e] for e in edges(s)))
-    vd = vertex_data(s)
-    ts = Dict{eltype(vertices(s)), ITensor}()
-    for v in vertices(s)
-        site_inds = isassigned(vd, v) ? vd[v] : Index[]
-        is = (site_inds..., (l[e] for e in incident_edges(s, v))...)
-        ts[v] = build_tensor(is)
-    end
-    return ITensorNetwork(ts)
-end
 
-# Build an ITensor network on a graph specified by the inds network `s`.
-# `link_space` sets the bond dimension. Entries are drawn from a standard
-# normal distribution.
+# --- random_tensornetwork ----------------------------------------------------
+
+# Core: at each vertex of `graph`, place an `itensor(randn(rng, eltype, ...), inds_v)`
+# whose inds are `siteinds[v]` concatenated with one fresh `Index(link_space, "Link")`
+# per incident edge, shared with the other endpoint. `siteinds` is anything indexable
+# by vertex (`keys(siteinds)` matches `vertices(graph)`); use `Index[]` per vertex for
+# no site inds.
 function random_tensornetwork(
-        rng::AbstractRNG, eltype::Type, s::IndsNetwork; kwargs...
-    )
-    return _random_tensornetwork(
-        s, is -> itensor(randn(rng, eltype, dim.(is)...), is...); kwargs...
+        rng::AbstractRNG, eltype::Type, graph::AbstractGraph, siteinds; link_space = 1
     )
-end
-
-function random_tensornetwork(eltype::Type, s::IndsNetwork; kwargs...)
-    return random_tensornetwork(Random.default_rng(), eltype, s; kwargs...)
-end
-
-function random_tensornetwork(rng::AbstractRNG, s::IndsNetwork; kwargs...)
-    return random_tensornetwork(rng, Float64, s; kwargs...)
-end
-
-function random_tensornetwork(s::IndsNetwork; kwargs...)
-    return random_tensornetwork(Random.default_rng(), s; kwargs...)
-end
-
-@traitfn function random_tensornetwork(
-        rng::AbstractRNG, eltype::Type, g::::IsUnderlyingGraph; kwargs...
+    g = NamedGraph(graph)
+    links = Dict(e => Index(link_space, "Link") for e in edges(g))
+    links = merge(links, Dict(reverse(e) => links[e] for e in edges(g)))
+    tensors = Dict(
+        map(collect(vertices(g))) do v
+            link_v = [links[e] for e in incident_edges(g, v)]
+            inds_v = [siteinds[v]; link_v]
+            return v => itensor(randn(rng, eltype, dim.(inds_v)...), inds_v)
+        end
     )
-    return random_tensornetwork(rng, eltype, IndsNetwork(g); kwargs...)
-end
-
-@traitfn function random_tensornetwork(eltype::Type, g::::IsUnderlyingGraph; kwargs...)
-    return random_tensornetwork(Random.default_rng(), eltype, g; kwargs...)
+    return ITensorNetwork(tensors, g)
 end
 
-@traitfn function random_tensornetwork(rng::AbstractRNG, g::::IsUnderlyingGraph; kwargs...)
-    return random_tensornetwork(rng, Float64, g; kwargs...)
-end
-
-@traitfn function random_tensornetwork(g::::IsUnderlyingGraph; kwargs...)
-    return random_tensornetwork(Random.default_rng(), g; kwargs...)
-end
-
-# Build an ITensor network on a graph specified by the inds network `s`.
-# Entries are drawn from the supplied `Distribution`.
+# `IndsNetwork`: extract site inds (`Index[]` where unassigned).
 function random_tensornetwork(
-        rng::AbstractRNG, distribution::Distribution, s::IndsNetwork; kwargs...
+        rng::AbstractRNG, eltype::Type, s::IndsNetwork; kwargs...
     )
-    return _random_tensornetwork(
-        s, is -> itensor(rand(rng, distribution, dim.(is)...), is...); kwargs...
+    siteinds = Dict(
+        v => isassigned(vertex_data(s), v) ? s[v] : Index[] for v in vertices(s)
     )
+    return random_tensornetwork(rng, eltype, underlying_graph(s), siteinds; kwargs...)
 end
 
-function random_tensornetwork(distribution::Distribution, s::IndsNetwork; kwargs...)
-    return random_tensornetwork(Random.default_rng(), distribution, s; kwargs...)
-end
-
-@traitfn function random_tensornetwork(
-        rng::AbstractRNG, distribution::Distribution, g::::IsUnderlyingGraph; kwargs...
+# Plain graph: no site inds.
+function random_tensornetwork(
+        rng::AbstractRNG, eltype::Type, g::AbstractGraph; kwargs...
     )
-    return random_tensornetwork(rng, distribution, IndsNetwork(g); kwargs...)
-end
-
-@traitfn function random_tensornetwork(
-        distribution::Distribution, g::::IsUnderlyingGraph; kwargs...
+    return random_tensornetwork(
+        rng,
+        eltype,
+        g,
+        Dict(v => Index[] for v in vertices(g));
+        kwargs...
     )
-    return random_tensornetwork(Random.default_rng(), distribution, g; kwargs...)
 end
 
-function random_ttn(args...; kwargs...)
-    return normalize(
-        ITensorNetworks._TreeTensorNetwork(random_tensornetwork(args...; kwargs...))
+# RNG / eltype / both defaults
+function random_tensornetwork(rng::AbstractRNG, sites; kwargs...)
+    return random_tensornetwork(rng, Float64, sites; kwargs...)
+end
+function random_tensornetwork(eltype::Type, sites; kwargs...)
+    return random_tensornetwork(Random.default_rng(), eltype, sites; kwargs...)
+end
+function random_tensornetwork(sites; kwargs...)
+    return random_tensornetwork(Random.default_rng(), Float64, sites; kwargs...)
+end
+
+# --- productstate -------------------------------------------------------------
+
+# Thread a fresh trivial-QN dim-1 link Index between `src(edge)` and `dst(edge)`
+# by outer-producting `onehot` basis vectors into each endpoint tensor. This is
+# the productstate-specific analog of `ITensorNetworks.add_edge!` — the latter
+# uses QR, which would push site-state QN flux into the link and leave BP's
+# `default_message` (single-index `delta`) with no compatible block. `onehot`
+# defaults to `Float64`, so we typecast it to `elt` to keep `productstate(elt,
+# ...)` element-type-preserving.
+function _add_edge!(elt::Type, tn, edge)
+    iₑ = Index(trivial_space(tn), "Link")
+    X = ITensors.convert_eltype(elt, onehot(iₑ => 1))
+    tn[src(edge)] *= X
+    tn[dst(edge)] *= dag(X)
+    return tn
+end
+
+# Build a product-state ITensorNetwork: start from a site-only TN (one tensor
+# per vertex with just the on-site state vector, looked up via
+# `ITensors.state(name, site_index)`), then add each edge from the original
+# IndsNetwork via `_add_edge!`. `state` is anything indexable by vertex (dict,
+# dictionary, array, ...); the `Function` method just materializes a Dict first.
+function productstate(elt::Type, state::Function, s::IndsNetwork)
+    return productstate(elt, Dict(v => state(v) for v in vertices(s)), s)
+end
+function productstate(elt::Type, state, s::IndsNetwork)
+    g = NamedGraph(collect(vertices(s)))
+    tensors = Dict(
+        map(collect(vertices(s))) do v
+            site_v = isassigned(vertex_data(s), v) ? s[v] : Index[]
+            t = ITensors.state(state[v], only(site_v))
+            return v => ITensors.convert_eltype(elt, t)
+        end
     )
+    tn = ITensorNetwork(tensors, g)
+    for e in edges(s)
+        _add_edge!(elt, tn, e)
+    end
+    return tn
 end
+productstate(state, s::IndsNetwork) = productstate(Float64, state, s)
 
-function random_mps(args...; kwargs...)
-    return random_ttn(args...; kwargs...)
-end
-
-function random_mps(f, is::Vector{<:Index}; kwargs...)
-    return random_mps(f, ITensorNetworks.path_indsnetwork(is); kwargs...)
-end
-
-function random_mps(s::Vector{<:Index}; kwargs...)
-    return random_mps(ITensorNetworks.path_indsnetwork(s); kwargs...)
-end
+# --- ModelHamiltonians --------------------------------------------------------
 
 # Small grab-bag of model-Hamiltonian builders used in regression tests. Kept
 # in a submodule so call sites remain `ModelHamiltonians.ising(g; h = ...)` etc.