diff --git a/Project.toml b/Project.toml
index 68ece85cd..fc4ad69c9 100644
--- a/Project.toml
+++ b/Project.toml
@@ -24,6 +24,7 @@ Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 TensorKit = "07d1fe3e-3e46-537d-9eac-e9e13d0d4cec"
+TensorKitTensors = "41b62e7d-e9d1-4e23-942c-79a97adf954b"
 TensorOperations = "6aa20fa7-93e2-5fca-9bc0-fbd0db3c71a2"
 TupleTools = "9d95972d-f1c8-5527-a6e0-b4b365fa01f6"
 VectorInterface = "409d34a3-91d5-4945-b6ec-7529ddf182d8"
@@ -47,6 +48,7 @@ Printf = "1"
 Random = "1"
 Statistics = "1"
 TensorKit = "0.16.5, 0.17"
+TensorKitTensors = "0.2.5"
 TensorOperations = "5"
 TupleTools = "1.6.0"
 VectorInterface = "0.4, 0.5, 0.6"
diff --git a/docs/Project.toml b/docs/Project.toml
index aecf8003d..e8d21f68b 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -7,6 +7,7 @@ MPSKit = "bb1c41ca-d63c-52ed-829e-0820dda26502"
 MPSKitModels = "ca635005-6f8c-4cd1-b51d-8491250ef2ab"
 PEPSKit = "52969e89-939e-4361-9b68-9bc7cde4bdeb"
 TensorKit = "07d1fe3e-3e46-537d-9eac-e9e13d0d4cec"
+TensorKitTensors = "41b62e7d-e9d1-4e23-942c-79a97adf954b"
 
 [compat]
 Documenter = "1"
diff --git a/docs/make.jl b/docs/make.jl
index 636330ec3..9f82c9980 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -24,6 +24,7 @@ links = InterLinks(
     "KrylovKit" => "https://jutho.github.io/KrylovKit.jl/stable/",
     "MPSKit" => "https://quantumkithub.github.io/MPSKit.jl/stable/",
     "MPSKitModels" => "https://quantumkithub.github.io/MPSKitModels.jl/dev/",
+    "TensorKitTensors" => "https://quantumkithub.github.io/TensorKitTensors.jl/stable/",
     # "Zygote" => "https://fluxml.ai/Zygote.jl/stable/",
     "ChainRulesCore" => "https://juliadiff.org/ChainRulesCore.jl/stable/",
     "MatrixAlgebraKit" => "https://quantumkithub.github.io/MatrixAlgebraKit.jl/stable/",
diff --git a/docs/src/examples/bose_hubbard/index.md b/docs/src/examples/bose_hubbard/index.md
index ad3dc7c1c..19a1681cf 100644
--- a/docs/src/examples/bose_hubbard/index.md
+++ b/docs/src/examples/bose_hubbard/index.md
@@ -27,8 +27,7 @@ Random.seed!(2928528935);
 ## Defining the model
 
 We will construct the Bose-Hubbard model Hamiltonian through the
-[`bose_hubbard_model`](https://quantumkithub.github.io/MPSKitModels.jl/dev/man/models/#MPSKitModels.bose_hubbard_model),
-function from MPSKitModels as reexported by PEPSKit. We'll simulate the model in its
+[`bose_hubbard_model`](@ref) defined in PEPSKit. We'll simulate the model in its
 Mott-insulating phase where the ratio $U/t$ is large, since in this phase we expect the
 ground state to be well approximated by a PEPS with a manifest global $U(1)$ symmetry.
 Furthermore, we'll impose a cutoff at 2 bosons per site, set the chemical potential to zero
diff --git a/docs/src/examples/bose_hubbard/main.ipynb b/docs/src/examples/bose_hubbard/main.ipynb
index 9c4ba4e3d..5f2e82a13 100644
--- a/docs/src/examples/bose_hubbard/main.ipynb
+++ b/docs/src/examples/bose_hubbard/main.ipynb
@@ -39,17 +39,7 @@
   {
    "cell_type": "markdown",
    "metadata": {},
-   "source": [
-    "## Defining the model\n",
-    "\n",
-    "We will construct the Bose-Hubbard model Hamiltonian through the\n",
-    "[`bose_hubbard_model`](https://quantumkithub.github.io/MPSKitModels.jl/dev/man/models/#MPSKitModels.bose_hubbard_model),\n",
-    "function from MPSKitModels as reexported by PEPSKit. We'll simulate the model in its\n",
-    "Mott-insulating phase where the ratio $U/t$ is large, since in this phase we expect the\n",
-    "ground state to be well approximated by a PEPS with a manifest global $U(1)$ symmetry.\n",
-    "Furthermore, we'll impose a cutoff at 2 bosons per site, set the chemical potential to zero\n",
-    "and use a simple $1 \\times 1$ unit cell:"
-   ]
+   "source": "## Defining the model\n\nWe will construct the Bose-Hubbard model Hamiltonian through the\n[`bose_hubbard_model`](@ref) defined in PEPSKit. We'll simulate the model in its\nMott-insulating phase where the ratio $U/t$ is large, since in this phase we expect the\nground state to be well approximated by a PEPS with a manifest global $U(1)$ symmetry.\nFurthermore, we'll impose a cutoff at 2 bosons per site, set the chemical potential to zero\nand use a simple $1 \\times 1$ unit cell:"
   },
   {
    "cell_type": "code",
diff --git a/docs/src/examples/c4v_ctmrg/index.md b/docs/src/examples/c4v_ctmrg/index.md
index 20a87b1be..fda303e52 100644
--- a/docs/src/examples/c4v_ctmrg/index.md
+++ b/docs/src/examples/c4v_ctmrg/index.md
@@ -40,7 +40,7 @@ by evaluating only half of the terms and multiplying by 2. In practice, we imple
 using a specialized [`LocalOperator`](@ref) that contains only the relevant terms:
 
 ````julia
-using MPSKitModels: S_xx, S_yy, S_zz
+import TensorKitTensors.SpinOperators as SO
 
 # Heisenberg model assuming C4v symmetric PEPS and environment, which only evaluates necessary term
 function heisenberg_XYZ_c4v(lattice::InfiniteSquare; kwargs...)
@@ -52,9 +52,9 @@ function heisenberg_XYZ_c4v(
     )
     @assert size(lattice) == (1, 1) "only trivial unit cells supported by C4v-symmetric Hamiltonians"
     term =
-        S_xx(T, S; spin = spin) * Jx +
-        S_yy(T, S; spin = spin) * Jy +
-        S_zz(T, S; spin = spin) * Jz
+        SO.S_x_S_x(T, S; spin = spin) * Jx +
+        SO.S_y_S_y(T, S; spin = spin) * Jy +
+        SO.S_z_S_z(T, S; spin = spin) * Jz
     spaces = fill(domain(term)[1], (1, 1))
     return LocalOperator( # horizontal and vertical contributions are identical
         spaces, (CartesianIndex(1, 1), CartesianIndex(1, 2)) => 2 * term
diff --git a/docs/src/examples/c4v_ctmrg/main.ipynb b/docs/src/examples/c4v_ctmrg/main.ipynb
index edc65df99..4cbc2f0c5 100644
--- a/docs/src/examples/c4v_ctmrg/main.ipynb
+++ b/docs/src/examples/c4v_ctmrg/main.ipynb
@@ -61,28 +61,7 @@
    "execution_count": null,
    "metadata": {},
    "outputs": [],
-   "source": [
-    "using MPSKitModels: S_xx, S_yy, S_zz\n",
-    "\n",
-    "# Heisenberg model assuming C4v symmetric PEPS and environment, which only evaluates necessary term\n",
-    "function heisenberg_XYZ_c4v(lattice::InfiniteSquare; kwargs...)\n",
-    "    return heisenberg_XYZ_c4v(ComplexF64, Trivial, lattice; kwargs...)\n",
-    "end\n",
-    "function heisenberg_XYZ_c4v(\n",
-    "        T::Type{<:Number}, S::Type{<:Sector}, lattice::InfiniteSquare;\n",
-    "        Jx = -1.0, Jy = 1.0, Jz = -1.0, spin = 1 // 2,\n",
-    "    )\n",
-    "    @assert size(lattice) == (1, 1) \"only trivial unit cells supported by C4v-symmetric Hamiltonians\"\n",
-    "    term =\n",
-    "        S_xx(T, S; spin = spin) * Jx +\n",
-    "        S_yy(T, S; spin = spin) * Jy +\n",
-    "        S_zz(T, S; spin = spin) * Jz\n",
-    "    spaces = fill(domain(term)[1], (1, 1))\n",
-    "    return LocalOperator( # horizontal and vertical contributions are identical\n",
-    "        spaces, (CartesianIndex(1, 1), CartesianIndex(1, 2)) => 2 * term\n",
-    "    )\n",
-    "end;"
-   ]
+   "source": "import TensorKitTensors.SpinOperators as SO\n\n# Heisenberg model assuming C4v symmetric PEPS and environment, which only evaluates necessary term\nfunction heisenberg_XYZ_c4v(lattice::InfiniteSquare; kwargs...)\n    return heisenberg_XYZ_c4v(ComplexF64, Trivial, lattice; kwargs...)\nend\nfunction heisenberg_XYZ_c4v(\n        T::Type{<:Number}, S::Type{<:Sector}, lattice::InfiniteSquare;\n        Jx = -1.0, Jy = 1.0, Jz = -1.0, spin = 1 // 2,\n    )\n    @assert size(lattice) == (1, 1) \"only trivial unit cells supported by C4v-symmetric Hamiltonians\"\n    term =\n        SO.S_x_S_x(T, S; spin = spin) * Jx +\n        SO.S_y_S_y(T, S; spin = spin) * Jy +\n        SO.S_z_S_z(T, S; spin = spin) * Jz\n    spaces = fill(domain(term)[1], (1, 1))\n    return LocalOperator( # horizontal and vertical contributions are identical\n        spaces, (CartesianIndex(1, 1), CartesianIndex(1, 2)) => 2 * term\n    )\nend;"
   },
   {
    "cell_type": "markdown",
diff --git a/docs/src/examples/heisenberg/index.md b/docs/src/examples/heisenberg/index.md
index 90607c846..24438b941 100644
--- a/docs/src/examples/heisenberg/index.md
+++ b/docs/src/examples/heisenberg/index.md
@@ -40,8 +40,7 @@ using TensorKit, PEPSKit
 ## Defining the Heisenberg Hamiltonian
 
 To create the sublattice rotated Heisenberg Hamiltonian on an infinite square lattice, we use
-the `heisenberg_XYZ` method from [MPSKitModels](https://quantumkithub.github.io/MPSKitModels.jl/dev/)
-which is redefined for the `InfiniteSquare` and reexported in PEPSKit:
+the `heisenberg_XYZ` model defined in PEPSKit for the `InfiniteSquare` lattice:
 
 ````julia
 H = heisenberg_XYZ(InfiniteSquare(); Jx = -1, Jy = 1, Jz = -1)
diff --git a/docs/src/examples/heisenberg/main.ipynb b/docs/src/examples/heisenberg/main.ipynb
index c8d0bd6f5..95a37a5f6 100644
--- a/docs/src/examples/heisenberg/main.ipynb
+++ b/docs/src/examples/heisenberg/main.ipynb
@@ -62,13 +62,7 @@
   {
    "cell_type": "markdown",
    "metadata": {},
-   "source": [
-    "## Defining the Heisenberg Hamiltonian\n",
-    "\n",
-    "To create the sublattice rotated Heisenberg Hamiltonian on an infinite square lattice, we use\n",
-    "the `heisenberg_XYZ` method from [MPSKitModels](https://quantumkithub.github.io/MPSKitModels.jl/dev/)\n",
-    "which is redefined for the `InfiniteSquare` and reexported in PEPSKit:"
-   ]
+   "source": "## Defining the Heisenberg Hamiltonian\n\nTo create the sublattice rotated Heisenberg Hamiltonian on an infinite square lattice, we use\nthe `heisenberg_XYZ` model defined in PEPSKit for the `InfiniteSquare` lattice:"
   },
   {
    "cell_type": "code",
diff --git a/docs/src/examples/heisenberg_su/index.md b/docs/src/examples/heisenberg_su/index.md
index 7f6056e56..056f4fa53 100644
--- a/docs/src/examples/heisenberg_su/index.md
+++ b/docs/src/examples/heisenberg_su/index.md
@@ -26,7 +26,7 @@ Let's get started by seeding the RNG and importing all required modules:
 using Random
 import Statistics: mean
 using TensorKit, PEPSKit
-import MPSKitModels: S_x, S_y, S_z, S_exchange
+import TensorKitTensors.SpinOperators as SO
 Random.seed!(0);
 ````
 
@@ -169,9 +169,9 @@ function compute_mags(peps::InfinitePEPS, env::CTMRGEnv)
     # detect symmetry on physical axis
     symm = sectortype(space(peps.A[1, 1]))
     if symm == Trivial
-        S_ops = real.([S_x(symm), im * S_y(symm), S_z(symm)])
+        S_ops = real.([SO.S_x(symm), im * SO.S_y(symm), SO.S_z(symm)])
     elseif symm == U1Irrep
-        S_ops = real.([S_z(symm)]) ## only Sz preserves <Sz>
+        S_ops = real.([SO.S_z(symm)]) ## only Sz preserves <Sz>
     end
 
     return map(Iterators.product(axes(peps, 1), axes(peps, 2), S_ops)) do (r, c, S)
diff --git a/docs/src/examples/heisenberg_su/main.ipynb b/docs/src/examples/heisenberg_su/main.ipynb
index 1a6c0ff39..a61c9ebf1 100644
--- a/docs/src/examples/heisenberg_su/main.ipynb
+++ b/docs/src/examples/heisenberg_su/main.ipynb
@@ -34,13 +34,7 @@
    "execution_count": null,
    "metadata": {},
    "outputs": [],
-   "source": [
-    "using Random\n",
-    "import Statistics: mean\n",
-    "using TensorKit, PEPSKit\n",
-    "import MPSKitModels: S_x, S_y, S_z, S_exchange\n",
-    "Random.seed!(0);"
-   ]
+   "source": "using Random\nimport Statistics: mean\nusing TensorKit, PEPSKit\nimport TensorKitTensors.SpinOperators as SO\nRandom.seed!(0);"
   },
   {
    "cell_type": "markdown",
@@ -175,30 +169,7 @@
    "execution_count": null,
    "metadata": {},
    "outputs": [],
-   "source": [
-    "function compute_mags(peps::InfinitePEPS, env::CTMRGEnv)\n",
-    "    lattice = collect(space(t, 1) for t in peps.A)\n",
-    "\n",
-    "    # detect symmetry on physical axis\n",
-    "    symm = sectortype(space(peps.A[1, 1]))\n",
-    "    if symm == Trivial\n",
-    "        S_ops = real.([S_x(symm), im * S_y(symm), S_z(symm)])\n",
-    "    elseif symm == U1Irrep\n",
-    "        S_ops = real.([S_z(symm)]) ## only Sz preserves <Sz>\n",
-    "    end\n",
-    "\n",
-    "    return map(Iterators.product(axes(peps, 1), axes(peps, 2), S_ops)) do (r, c, S)\n",
-    "        expectation_value(peps, LocalOperator(lattice, (CartesianIndex(r, c),) => S), env)\n",
-    "    end\n",
-    "end\n",
-    "\n",
-    "E = expectation_value(peps, H, env) / (Nr * Nc)\n",
-    "Ms = compute_mags(peps, env)\n",
-    "M_norms = map(\n",
-    "    rc -> norm(Ms[rc[1], rc[2], :]), Iterators.product(axes(peps, 1), axes(peps, 2))\n",
-    ")\n",
-    "@show E Ms M_norms;"
-   ]
+   "source": "function compute_mags(peps::InfinitePEPS, env::CTMRGEnv)\n    lattice = collect(space(t, 1) for t in peps.A)\n\n    # detect symmetry on physical axis\n    symm = sectortype(space(peps.A[1, 1]))\n    if symm == Trivial\n        S_ops = real.([SO.S_x(symm), im * SO.S_y(symm), SO.S_z(symm)])\n    elseif symm == U1Irrep\n        S_ops = real.([SO.S_z(symm)]) ## only Sz preserves <Sz>\n    end\n\n    return map(Iterators.product(axes(peps, 1), axes(peps, 2), S_ops)) do (r, c, S)\n        expectation_value(peps, LocalOperator(lattice, (CartesianIndex(r, c),) => S), env)\n    end\nend\n\nE = expectation_value(peps, H, env) / (Nr * Nc)\nMs = compute_mags(peps, env)\nM_norms = map(\n    rc -> norm(Ms[rc[1], rc[2], :]), Iterators.product(axes(peps, 1), axes(peps, 2))\n)\n@show E Ms M_norms;"
   },
   {
    "cell_type": "markdown",
diff --git a/docs/src/man/models.md b/docs/src/man/models.md
index 5539fb6f8..ca9666435 100644
--- a/docs/src/man/models.md
+++ b/docs/src/man/models.md
@@ -1,29 +1,32 @@
 # Models
 
-PEPSKit implements physical models through the [MPSKitModels.jl](https://quantumkithub.github.io/MPSKitModels.jl/dev/) package as [`PEPSKit.LocalOperator`](@ref) structs.
+PEPSKit implements physical models as [`PEPSKit.LocalOperator`](@ref) structs.
 Here, we want to explain how users can define their own Hamiltonians and provide a list of
 already implemented models.
 
 ## Implementing custom models
 
-In order to define custom Hamiltonians, we leverage several of the useful tools provided in MPSKitModels.
-In particular, we use many of the pre-defined [operators](@extref MPSKitModels Operators), which is especially useful when defining models with symmetric and fermionic tensors, since most of these operators can take a symmetry as an argument, returning the appropriate symmetric `TensorMap`.
+In order to define custom Hamiltonians, we leverage the operator building blocks provided by
+[TensorKitTensors.jl](https://quantumkithub.github.io/TensorKitTensors.jl/stable/), which offers
+pre-defined symmetric tensors for spin, boson, fermion, and Hubbard systems.
 In order to specify the lattice on which the Hamiltonian is defined, we construct two-dimensional lattices as subtypes of [`MPSKitModels.AbstractLattice`](@extref).
 Note that so far, all models are defined on infinite square lattices, see [`InfiniteSquare`](@ref), but in the future, we plan to support other lattice geometries as well.
 In order to specify tensors acting on particular lattice sites, there are a couple of handy methods that we want to point to: see `vertices`, `nearest_neighbors` and `next_nearest_neighbors` defined [here](https://github.com/QuantumKitHub/PEPSKit.jl/blob/master/src/operators/lattices/squarelattice.jl).
 
-For a simple example on how to implement a custom model, let's look at the implementation of the [`MPSKitModels.transverse_field_ising`](@extref) model:
+For a simple example on how to implement a custom model, let's look at the implementation of the
+[`transverse_field_ising`](@ref) model:
 
 ```julia
-function MPSKitModels.transverse_field_ising(
+import TensorKitTensors.SpinOperators as SO
+function transverse_field_ising(
     T::Type{<:Number},
     S::Union{Type{Trivial},Type{Z2Irrep}},
     lattice::InfiniteSquare;
     J=1.0,
     g=1.0,
 )
-    ZZ = rmul!(σᶻᶻ(T, S), -J)
-    X = rmul!(σˣ(T, S), g * -J)
+    ZZ = rmul!(SO.S_z_S_z(T, S), -4 * J)
+    X = rmul!(SO.σˣ(T, S), g * -J)
     spaces = fill(domain(X)[1], (lattice.Nrows, lattice.Ncols))
     return LocalOperator(
         spaces,
@@ -35,7 +38,7 @@ end
 
 This provides a good recipe for defining a model:
 
-1. Define the locally-acting tensors as `TensorMap`s.
+1. Define the locally-acting tensors as `TensorMap`s using `TensorKitTensors` operators.
 2. Construct a matrix of the physical spaces these `TensorMap`s act on based on the lattice geometry.
 3. Return a `LocalOperator` where we specify on which sites (e.g. on-site, nearest neighbor, etc.) the local tensors act.
 
@@ -43,20 +46,21 @@ For more model implementations, check the [PEPSKit repository](https://github.co
 
 ## Implemented models
 
-While PEPSKit provides an interface for specifying custom Hamiltonians, it also provides a number of [pre-defined models](https://github.com/QuantumKitHub/PEPSKit.jl/blob/master/src/operators/models.jl). Some of these are models already defined in [MPSKitModels](@extref MPSKitModels Models), which are overloaded for two-dimensional lattices and re-exported, but there are new additions as well. The following models are provided:
+PEPSKit provides a number of pre-defined models. The following model constructors are available
+and can be used directly with an [`InfiniteSquare`](@ref) lattice:
 
-### MPSKitModels.jl models
+### Models inherited from MPSKitModels.jl
 
-```@docs
-MPSKitModels.transverse_field_ising
-MPSKitModels.heisenberg_XYZ
-MPSKitModels.heisenberg_XXZ
-MPSKitModels.hubbard_model
-MPSKitModels.bose_hubbard_model
-MPSKitModels.tj_model
+```@docs; canonical=false
+transverse_field_ising(::InfiniteSquare)
+heisenberg_XYZ(::InfiniteSquare)
+heisenberg_XXZ(::InfiniteSquare)
+hubbard_model(::InfiniteSquare)
+bose_hubbard_model(::InfiniteSquare)
+tj_model(::InfiniteSquare)
 ```
 
-### PEPSKit.jl models
+### Models introduced by PEPSKit.jl
 
 ```@docs
 j1_j2_model
diff --git a/docs/src/man/precompilation.md b/docs/src/man/precompilation.md
index 2265c3603..85c1fe2e6 100644
--- a/docs/src/man/precompilation.md
+++ b/docs/src/man/precompilation.md
@@ -40,7 +40,7 @@ module Startup
 using Random
 using TensorKit, KrylovKit, OptimKit
 using ChainRulesCore, Zygote
-using MPSKit, MPSKitModels
+using MPSKit
 using PEPSKit
 using PrecompileTools
 
diff --git a/examples/Project.toml b/examples/Project.toml
index b66b36c4a..92bfaf4a1 100644
--- a/examples/Project.toml
+++ b/examples/Project.toml
@@ -3,11 +3,11 @@ KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
 MPSKit = "bb1c41ca-d63c-52ed-829e-0820dda26502"
-MPSKitModels = "ca635005-6f8c-4cd1-b51d-8491250ef2ab"
 OptimKit = "77e91f04-9b3b-57a6-a776-40b61faaebe0"
 PEPSKit = "52969e89-939e-4361-9b68-9bc7cde4bdeb"
 QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 TensorKit = "07d1fe3e-3e46-537d-9eac-e9e13d0d4cec"
+TensorKitTensors = "41b62e7d-e9d1-4e23-942c-79a97adf954b"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
diff --git a/examples/bose_hubbard/main.jl b/examples/bose_hubbard/main.jl
index 05af563a6..813a29194 100644
--- a/examples/bose_hubbard/main.jl
+++ b/examples/bose_hubbard/main.jl
@@ -20,8 +20,7 @@ md"""
 ## Defining the model
 
 We will construct the Bose-Hubbard model Hamiltonian through the
-[`bose_hubbard_model`](https://quantumkithub.github.io/MPSKitModels.jl/dev/man/models/#MPSKitModels.bose_hubbard_model),
-function from MPSKitModels as reexported by PEPSKit. We'll simulate the model in its
+[`bose_hubbard_model`](@ref) defined in PEPSKit. We'll simulate the model in its
 Mott-insulating phase where the ratio $U/t$ is large, since in this phase we expect the
 ground state to be well approximated by a PEPS with a manifest global $U(1)$ symmetry.
 Furthermore, we'll impose a cutoff at 2 bosons per site, set the chemical potential to zero
diff --git a/examples/c4v_ctmrg/main.jl b/examples/c4v_ctmrg/main.jl
index ec2d47703..6558494ea 100644
--- a/examples/c4v_ctmrg/main.jl
+++ b/examples/c4v_ctmrg/main.jl
@@ -33,7 +33,7 @@ by evaluating only half of the terms and multiplying by 2. In practice, we imple
 using a specialized [`LocalOperator`](@ref) that contains only the relevant terms:
 """
 
-using MPSKitModels: S_xx, S_yy, S_zz
+import TensorKitTensors.SpinOperators as SO
 
 ## Heisenberg model assuming C4v symmetric PEPS and environment, which only evaluates necessary term
 function heisenberg_XYZ_c4v(lattice::InfiniteSquare; kwargs...)
@@ -45,9 +45,9 @@ function heisenberg_XYZ_c4v(
     )
     @assert size(lattice) == (1, 1) "only trivial unit cells supported by C4v-symmetric Hamiltonians"
     term =
-        S_xx(T, S; spin = spin) * Jx +
-        S_yy(T, S; spin = spin) * Jy +
-        S_zz(T, S; spin = spin) * Jz
+        SO.S_x_S_x(T, S; spin = spin) * Jx +
+        SO.S_y_S_y(T, S; spin = spin) * Jy +
+        SO.S_z_S_z(T, S; spin = spin) * Jz
     spaces = fill(domain(term)[1], (1, 1))
     return LocalOperator( # horizontal and vertical contributions are identical
         spaces, (CartesianIndex(1, 1), CartesianIndex(1, 2)) => 2 * term
diff --git a/examples/heisenberg/main.jl b/examples/heisenberg/main.jl
index 5f571d4ec..948a8588a 100644
--- a/examples/heisenberg/main.jl
+++ b/examples/heisenberg/main.jl
@@ -33,8 +33,7 @@ md"""
 ## Defining the Heisenberg Hamiltonian
 
 To create the sublattice rotated Heisenberg Hamiltonian on an infinite square lattice, we use
-the `heisenberg_XYZ` method from [MPSKitModels](https://quantumkithub.github.io/MPSKitModels.jl/dev/)
-which is redefined for the `InfiniteSquare` and reexported in PEPSKit:
+the `heisenberg_XYZ` model defined in PEPSKit for the `InfiniteSquare` lattice:
 """
 
 H = heisenberg_XYZ(InfiniteSquare(); Jx = -1, Jy = 1, Jz = -1)
diff --git a/examples/heisenberg_su/main.jl b/examples/heisenberg_su/main.jl
index 00d0eb100..7beec5e3c 100644
--- a/examples/heisenberg_su/main.jl
+++ b/examples/heisenberg_su/main.jl
@@ -19,7 +19,7 @@ Let's get started by seeding the RNG and importing all required modules:
 using Random
 import Statistics: mean
 using TensorKit, PEPSKit
-import MPSKitModels: S_x, S_y, S_z, S_exchange
+import TensorKitTensors.SpinOperators as SO
 Random.seed!(0);
 
 md"""
@@ -112,9 +112,9 @@ function compute_mags(peps::InfinitePEPS, env::CTMRGEnv)
     ## detect symmetry on physical axis
     symm = sectortype(space(peps.A[1, 1]))
     if symm == Trivial
-        S_ops = real.([S_x(symm), im * S_y(symm), S_z(symm)])
+        S_ops = real.([SO.S_x(symm), im * SO.S_y(symm), SO.S_z(symm)])
     elseif symm == U1Irrep
-        S_ops = real.([S_z(symm)]) ## only Sz preserves <Sz>
+        S_ops = real.([SO.S_z(symm)]) ## only Sz preserves <Sz>
     end
 
     return map(Iterators.product(axes(peps, 1), axes(peps, 2), S_ops)) do (r, c, S)
diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl
index a6e199c54..7a1ba2677 100644
--- a/src/PEPSKit.jl
+++ b/src/PEPSKit.jl
@@ -33,7 +33,15 @@ using MPSKit: InfiniteEnvironments
 import MPSKit: tensorexpr, leading_boundary, loginit!, logiter!, logfinish!, logcancel!, physicalspace
 import MPSKit: infinite_temperature_density_matrix
 
+import TensorKitTensors.SpinOperators as SO
+import TensorKitTensors.FermionOperators as FO
+import TensorKitTensors.HubbardOperators as HO
+import TensorKitTensors.TJOperators as TJ
+import TensorKitTensors.BosonOperators as BO
 using MPSKitModels
+import MPSKitModels:
+    transverse_field_ising, heisenberg_XYZ, heisenberg_XXZ,
+    hubbard_model, bose_hubbard_model, tj_model
 using FiniteDifferences
 using OhMyThreads: tmap, tmap!
 using DocStringExtensions
diff --git a/src/operators/models.jl b/src/operators/models.jl
index 6e1d927b9..0a9c6e051 100644
--- a/src/operators/models.jl
+++ b/src/operators/models.jl
@@ -22,15 +22,33 @@ function nearest_neighbour_hamiltonian(
 end
 
 #
-## MPSKitModels.jl re-exports
+## Model definitions
 #
 
-function MPSKitModels.transverse_field_ising(
+"""
+    transverse_field_ising(
+        [T::Type{<:Number}], [S::Union{Type{Trivial}, Type{Z2Irrep}}],
+        lattice::InfiniteSquare; J=1.0, g=1.0
+    )
+
+`LocalOperator` for the [transverse-field Ising model](https://en.wikipedia.org/wiki/Transverse-field_Ising_model)
+Hamiltonian on the square lattice,
+```math
+H = -J\\left(\\sum_{\\langle i,j \\rangle} \\sigma^z_i \\sigma^z_j + g \\sum_{i} \\sigma^x_i \\right)
+```
+where ``\\sigma^i`` are the spin-1/2 Pauli operators.
+
+By default, it is defined with `Trivial` symmetry and `ComplexF64` entries.
+"""
+function transverse_field_ising(lattice::InfiniteSquare; kwargs...)
+    return transverse_field_ising(ComplexF64, Trivial, lattice; kwargs...)
+end
+function transverse_field_ising(
         T::Type{<:Number}, S::Union{Type{Trivial}, Type{Z2Irrep}}, lattice::InfiniteSquare;
         J = 1.0, g = 1.0,
     )
-    ZZ = rmul!(σᶻᶻ(T, S), -J)
-    X = rmul!(σˣ(T, S), g * -J)
+    ZZ = rmul!(SO.S_z_S_z(T, S), -4 * J)
+    X = rmul!(SO.σˣ(T, S), g * -J)
     spaces = fill(domain(X)[1], (lattice.Nrows, lattice.Ncols))
     return LocalOperator(
         spaces,
@@ -39,31 +57,60 @@ function MPSKitModels.transverse_field_ising(
     )
 end
 
-function MPSKitModels.heisenberg_XYZ(lattice::InfiniteSquare; kwargs...)
+"""
+    heisenberg_XYZ([T::Type{<:Number}], [S::Type{<:Sector}], lattice::InfiniteSquare;
+                   Jx=-1.0, Jy=1.0, Jz=-1.0, spin=1//2)
+
+`LocalOperator` for the XYZ Heisenberg model Hamiltonian on the square lattice,
+```math
+H = \\sum_{\\langle i,j \\rangle} \\left( J_x S_i^x S_j^x + J_y S_i^y S_j^y + J_z S_i^z S_j^z \\right)
+```
+
+By default, it is defined with `Trivial` symmetry and `ComplexF64` entries.
+
+See also [`heisenberg_XXZ`](@ref).
+"""
+function heisenberg_XYZ(lattice::InfiniteSquare; kwargs...)
     return heisenberg_XYZ(ComplexF64, Trivial, lattice; kwargs...)
 end
-function MPSKitModels.heisenberg_XYZ(
+function heisenberg_XYZ(
         T::Type{<:Number}, S::Type{<:Sector}, lattice::InfiniteSquare;
         Jx = -1.0, Jy = 1.0, Jz = -1.0, spin = 1 // 2,
     )
     term =
-        rmul!(S_xx(T, S; spin = spin), Jx) +
-        rmul!(S_yy(T, S; spin = spin), Jy) +
-        rmul!(S_zz(T, S; spin = spin), Jz)
+        rmul!(SO.S_x_S_x(T, S; spin = spin), Jx) +
+        rmul!(SO.S_y_S_y(T, S; spin = spin), Jy) +
+        rmul!(SO.S_z_S_z(T, S; spin = spin), Jz)
     spaces = fill(domain(term)[1], (lattice.Nrows, lattice.Ncols))
     return LocalOperator(
         spaces, (neighbor => term for neighbor in nearest_neighbours(lattice))...
     )
 end
 
-function MPSKitModels.heisenberg_XXZ(
+"""
+    heisenberg_XXZ([T::Type{<:Number}], [S::Type{<:Sector}], lattice::InfiniteSquare;
+                   J=1.0, Delta=1.0, spin=1)
+
+`LocalOperator` for the XXZ Heisenberg model Hamiltonian on the square lattice,
+```math
+H = J \\sum_{\\langle i,j \\rangle} \\left( S_i^x S_j^x + S_i^y S_j^y + \\Delta S_i^z S_j^z \\right)
+```
+
+By default, it is defined with `Trivial` symmetry and `ComplexF64` entries.
+
+See also [`heisenberg_XYZ`](@ref).
+"""
+function heisenberg_XXZ(lattice::InfiniteSquare; kwargs...)
+    return heisenberg_XXZ(ComplexF64, Trivial, lattice; kwargs...)
+end
+function heisenberg_XXZ(
         T::Type{<:Number}, S::Type{<:Sector}, lattice::InfiniteSquare;
         J = 1.0, Delta = 1.0, spin = 1
     )
     h =
         J * (
-        (S_plusmin(T, S; spin = spin) + S_minplus(T, S; spin = spin)) / 2 +
-            Delta * S_zz(T, S; spin = spin)
+        (SO.S_plus_S_min(T, S; spin = spin) + SO.S_min_S_plus(T, S; spin = spin)) / 2 +
+            Delta * SO.S_z_S_z(T, S; spin = spin)
     )
     spaces = fill(domain(h)[1], (lattice.Nrows, lattice.Ncols))
     return LocalOperator(
@@ -71,33 +118,68 @@ function MPSKitModels.heisenberg_XXZ(
     )
 end
 
-function MPSKitModels.hubbard_model(
+"""
+    hubbard_model([T::Type{<:Number}], [particle_symmetry::Type{<:Sector}],
+                  [spin_symmetry::Type{<:Sector}], lattice::InfiniteSquare;
+                  t=1.0, U=1.0, mu=0.0, n=0)
+
+`LocalOperator` for the Fermi-Hubbard model Hamiltonian on the square lattice,
+```math
+H = -t \\sum_{\\langle i,j \\rangle} \\sum_{\\sigma}
+    \\left( e_{i,\\sigma}^\\dagger e_{j,\\sigma} + \\text{h.c.} \\right)
+    + U \\sum_i n_{i,\\uparrow} n_{i,\\downarrow} - \\mu \\sum_i n_i
+```
+where ``\\sigma \\in \\{\\uparrow, \\downarrow\\}`` is a spin index and ``n`` is the
+fermionic number operator.
+
+By default, it is defined without any symmetries and with `ComplexF64` entries.
+"""
+function hubbard_model(lattice::InfiniteSquare; kwargs...)
+    return hubbard_model(ComplexF64, Trivial, Trivial, lattice; kwargs...)
+end
+function hubbard_model(
         T::Type{<:Number}, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector},
         lattice::InfiniteSquare;
         t = 1.0, U = 1.0, mu = 0.0, n::Integer = 0,
     )
     # TODO: just add this
     @assert n == 0 "Currently no support for imposing a fixed particle number"
-    N = MPSKitModels.e_number(T, particle_symmetry, spin_symmetry)
+    N = HO.e_num(T, particle_symmetry, spin_symmetry)
     pspace = space(N, 1)
     unit = TensorKit.id(pspace)
-    hopping =
-        MPSKitModels.e⁺e⁻(T, particle_symmetry, spin_symmetry) +
-        MPSKitModels.e⁻e⁺(T, particle_symmetry, spin_symmetry)
-    interaction_term = MPSKitModels.nꜛnꜜ(T, particle_symmetry, spin_symmetry)
+    hopping = HO.e_hopping(T, particle_symmetry, spin_symmetry)
+    interaction_term = HO.ud_num(T, particle_symmetry, spin_symmetry)
     site_term = U * interaction_term - mu * N
     h = (-t) * hopping + (1 / 4) * (site_term ⊗ unit + unit ⊗ site_term)
     return nearest_neighbour_hamiltonian(fill(pspace, size(lattice)), h)
 end
 
-function MPSKitModels.bose_hubbard_model(
-        elt::Type{<:Number}, symmetry::Type{<:Sector}, lattice::InfiniteSquare;
+"""
+    bose_hubbard_model([elt::Type{<:Number}], [symmetry::Type{<:Sector}],
+                       lattice::InfiniteSquare; cutoff=5, t=1.0, U=1.0, mu=0.0, n=0)
+
+`LocalOperator` for the Bose-Hubbard model Hamiltonian on the square lattice,
+```math
+H = -t \\sum_{\\langle i,j \\rangle} \\left( a_i^\\dagger a_j + \\text{h.c.} \\right)
+    - \\mu \\sum_i N_i + \\frac{U}{2} \\sum_i N_i(N_i - 1)
+```
+where ``N = a^\\dagger a`` is the bosonic number operator.
+
+The Hilbert space is truncated such that at maximum `cutoff` bosons can occupy a single site.
+If `symmetry` is `U1Irrep`, a fixed (half-integer) particle number density `n` can be imposed.
+
+By default, it is defined with `Trivial` symmetry and `ComplexF64` entries.
+"""
+function bose_hubbard_model(lattice::InfiniteSquare; kwargs...)
+    return bose_hubbard_model(ComplexF64, Trivial, lattice; kwargs...)
+end
+function bose_hubbard_model(
+        T::Type{<:Number}, symmetry::Type{<:Sector}, lattice::InfiniteSquare;
         cutoff::Integer = 5, t = 1.0, U = 1.0, mu = 0.0, n::Integer = 0,
     )
-    hopping_term =
-        a_plusmin(elt, symmetry; cutoff = cutoff) + a_minplus(elt, symmetry; cutoff = cutoff)
-    N = a_number(elt, symmetry; cutoff = cutoff)
-    interaction_term = MPSKitModels.contract_onesite(N, N - id(domain(N)))
+    hopping_term = BO.b_hopping(T, symmetry; cutoff)
+    N = BO.b_num(T, symmetry; cutoff = cutoff)
+    interaction_term = N * (N - id(domain(N)))
 
     spaces = fill(space(N, 1), (lattice.Nrows, lattice.Ncols))
 
@@ -120,45 +202,52 @@ function MPSKitModels.bose_hubbard_model(
     return H
 end
 
-function MPSKitModels.tj_model(
+"""
+    tj_model([T::Type{<:Number}], [particle_symmetry::Type{<:Sector}],
+             [spin_symmetry::Type{<:Sector}], lattice::InfiniteSquare;
+             t=2.5, J=1.0, mu=0.0, slave_fermion=false)
+
+`LocalOperator` for the t-J model Hamiltonian on the square lattice,
+```math
+H = -t \\sum_{\\langle i,j \\rangle, \\sigma}
+    (\\tilde{e}^\\dagger_{i,\\sigma} \\tilde{e}_{j,\\sigma} + \\text{h.c.})
+    + J \\sum_{\\langle i,j \\rangle}(\\mathbf{S}_i \\cdot \\mathbf{S}_j - \\frac{1}{4} n_i n_j)
+    - \\mu \\sum_i n_i
+```
+where ``\\tilde{e}_{i,\\sigma}`` is the electron operator with spin ``\\sigma`` projected to
+the no-double-occupancy subspace.
+
+By default, it is defined without any symmetries and with `ComplexF64` entries.
+"""
+function tj_model(lattice::InfiniteSquare; kwargs...)
+    return tj_model(ComplexF64, Trivial, Trivial, lattice; kwargs...)
+end
+function tj_model(
         T::Type{<:Number}, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector},
         lattice::InfiniteSquare;
         t = 2.5, J = 1.0, mu = 0.0, slave_fermion::Bool = false,
     )
-    hopping =
-        TJOperators.e_plusmin(particle_symmetry, spin_symmetry; slave_fermion) +
-        TJOperators.e_minplus(particle_symmetry, spin_symmetry; slave_fermion)
-    num = TJOperators.e_number(particle_symmetry, spin_symmetry; slave_fermion)
-    heis =
-        TJOperators.S_exchange(particle_symmetry, spin_symmetry; slave_fermion) -
+    hopping = TJ.e_hopping(T, particle_symmetry, spin_symmetry; slave_fermion)
+    num = TJ.e_num(T, particle_symmetry, spin_symmetry; slave_fermion)
+    heis = TJ.S_exchange(T, particle_symmetry, spin_symmetry; slave_fermion) -
         (1 / 4) * (num ⊗ num)
     pspace = space(num, 1)
     unit = TensorKit.id(pspace)
     h = (-t) * hopping + J * heis - (mu / 4) * (num ⊗ unit + unit ⊗ num)
-    if T <: Real
-        h = real(h)
-    end
     return nearest_neighbour_hamiltonian(fill(pspace, size(lattice)), h)
 end
 
-#
-## Models not yet defined in MPSKitModels
-#
-
 """
-    j1_j2_model([elt::Type{T}, symm::Type{S},] lattice::InfiniteSquare;
+    j1_j2_model([T::Type{<:Number}, S::Type{<:Sector},] lattice::InfiniteSquare;
                 J1=1.0, J2=1.0, spin=1//2, sublattice=true)
 
-Square lattice ``J_1\\text{-}J_2`` model, defined by the Hamiltonian
-
+`LocalOperator` for the ``J_1\\text{-}J_2`` model Hamiltonian on the square lattice,
 ```math
 H = J_1 \\sum_{\\langle i,j \\rangle} \\vec{S}_i \\cdot \\vec{S}_j
 + J_2 \\sum_{\\langle\\langle i,j \\rangle\\rangle} \\vec{S}_i \\cdot \\vec{S}_j,
 ```
-
-where ``\\vec{S}_i = (S_i^x, S_i^y, S_i^z)``. We denote the nearest and next-nearest neighbor
-terms using ``\\langle i,j \\rangle`` and ``\\langle\\langle i,j \\rangle\\rangle``,
-respectively. The `sublattice` kwarg enables a single-site unit cell ground state via a
+where ``\\vec{S}_i = (S_i^x, S_i^y, S_i^z)``.
+The `sublattice` kwarg enables a single-site unit cell ground state via a
 unitary sublattice rotation.
 """
 function j1_j2_model(lattice::InfiniteSquare; kwargs...)
@@ -168,9 +257,9 @@ function j1_j2_model(
         T::Type{<:Number}, S::Type{<:Sector}, lattice::InfiniteSquare;
         J1 = 1.0, J2 = 1.0, spin = 1 // 2, sublattice = true,
     )
-    term_AA = S_exchange(T, S; spin)
+    term_AA = SO.S_exchange(T, S; spin)
     term_AB = if sublattice
-        -S_xx(T, S; spin) + S_yy(T, S; spin) - S_zz(T, S; spin)  # Apply sublattice rotation
+        -SO.S_x_S_x(T, S; spin) + SO.S_y_S_y(T, S; spin) - SO.S_z_S_z(T, S; spin)  # Apply sublattice rotation
     else
         term_AA
     end
@@ -185,15 +274,15 @@ end
 """
     pwave_superconductor([T=ComplexF64,] lattice::InfiniteSquare; t=1, μ=2, Δ=1)
 
-Square lattice ``p``-wave superconductor model, defined by the Hamiltonian
-
+`LocalOperator` for the ``(p - ip)``-wave superconductor Hamiltonian on the square lattice
 ```math
     H = -\\sum_{\\langle i,j \\rangle} \\left( t c_i^\\dagger c_j +
-    \\Delta c_i c_j + \\text{h.c.} \\right) - \\mu \\sum_i n_i,
+    \\Delta_{ij} c_i c_j + \\text{h.c.} \\right) - \\mu \\sum_i n_i,
 ```
-
-where ``t`` is the hopping amplitude, ``\\Delta`` specifies the superconducting gap, ``\\mu``
+where ``t`` is the hopping amplitude, ``\\Delta_{ij}`` specifies the superconducting gap, ``\\mu``
 is the chemical potential, and ``n_i = c_i^\\dagger c_i`` is the fermionic number operator.
+For ``(p - ip)``-wave, ``\\Delta_{ij} = \\Delta`` on horizontal bonds,
+and ``-i \\Delta`` on vertical bonds.
 """
 function pwave_superconductor(lattice::InfiniteSquare; kwargs...)
     return pwave_superconductor(ComplexF64, lattice; kwargs...)
@@ -202,22 +291,19 @@ function pwave_superconductor(
         T::Type{<:Number}, lattice::InfiniteSquare;
         t::Number = 1, μ::Number = 2, Δ::Number = 1
     )
-    physical_space = Vect[FermionParity](0 => 1, 1 => 1)
+    physical_space = FO.fermion_space(Trivial)
     spaces = fill(physical_space, (lattice.Nrows, lattice.Ncols))
+    # hopping and pairing operators
+    hopp = -t * FO.f_hopping(T, Trivial)
+    pair = FO.f_min_f_min(T, Trivial)
+    pair_x, pair_y = Δ * pair, -im * Δ * pair
 
     # on-site
-    h0 = zeros(T, physical_space ← physical_space)
-    block(h0, FermionParity(1)) .= -μ
-
+    h0 = -μ * FO.f_num(T, Trivial)
     # two-site (x-direction)
-    hx = zeros(T, physical_space^2 ← physical_space^2)
-    block(hx, FermionParity(0)) .= [0 -Δ; -Δ 0]
-    block(hx, FermionParity(1)) .= [0 -t; -t 0]
-
+    hx = hopp + (pair_x + pair_x')
     # two-site (y-direction)
-    hy = zeros(T, physical_space^2 ← physical_space^2)
-    block(hy, FermionParity(0)) .= [0 Δ * im; -Δ * im 0]
-    block(hy, FermionParity(1)) .= [0 -t; -t 0]
+    hy = hopp + (pair_y + pair_y')
 
     x_neighbors = filter(n -> n[2].I[2] > n[1].I[2], nearest_neighbours(lattice))
     y_neighbors = filter(n -> n[2].I[1] > n[1].I[1], nearest_neighbours(lattice))
diff --git a/test/Project.toml b/test/Project.toml
index 2742a2812..1961b9868 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -15,6 +15,7 @@ ParallelTestRunner = "d3525ed8-44d0-4b2c-a655-542cee43accc"
 QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 TensorKit = "07d1fe3e-3e46-537d-9eac-e9e13d0d4cec"
+TensorKitTensors = "41b62e7d-e9d1-4e23-942c-79a97adf954b"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TestExtras = "5ed8adda-3752-4e41-b88a-e8b09835ee3a"
 VectorInterface = "409d34a3-91d5-4945-b6ec-7529ddf182d8"
diff --git a/test/examples/heisenberg.jl b/test/examples/heisenberg.jl
index 2968e09c1..02361148a 100644
--- a/test/examples/heisenberg.jl
+++ b/test/examples/heisenberg.jl
@@ -6,7 +6,7 @@ using TensorKit
 using KrylovKit
 using OptimKit
 using PEPSKit: peps_normalize
-using MPSKitModels: S_xx, S_yy, S_zz
+import TensorKitTensors.SpinOperators as SO
 
 # initialize parameters
 Dbond = 2
@@ -26,9 +26,9 @@ function heisenberg_XYZ_c4v(
     )
     @assert size(lattice) == (1, 1) "only trivial unit cells supported by C4v-symmetric Hamiltonians"
     term =
-        rmul!(S_xx(T, S; spin = spin), Jx) +
-        rmul!(S_yy(T, S; spin = spin), Jy) +
-        rmul!(S_zz(T, S; spin = spin), Jz)
+        rmul!(SO.S_x_S_x(T, S; spin = spin), Jx) +
+        rmul!(SO.S_y_S_y(T, S; spin = spin), Jy) +
+        rmul!(SO.S_z_S_z(T, S; spin = spin), Jz)
     spaces = fill(domain(term)[1], (1, 1))
     return LocalOperator( # horizontal and vertical contributions are identical
         spaces, [CartesianIndex(1, 1), CartesianIndex(1, 2)] => 2 * term
diff --git a/test/timeevol/cluster_projectors.jl b/test/timeevol/cluster_projectors.jl
index f95b9043c..d1e5d27e0 100644
--- a/test/timeevol/cluster_projectors.jl
+++ b/test/timeevol/cluster_projectors.jl
@@ -3,7 +3,7 @@ using TensorKit
 using PEPSKit
 using LinearAlgebra
 using Random
-import MPSKitModels: hubbard_space
+import TensorKitTensors.HubbardOperators as HO
 using PEPSKit: sdiag_pow, _cluster_truncate!, _flip_virtuals!
 using MPSKit: GenericMPSTensor, MPSBondTensor
 
@@ -223,7 +223,7 @@ end
     ctmrg_tol = 1.0e-9
     Random.seed!(1459)
     # with U(1) spin rotation symmetry
-    Pspace = hubbard_space(Trivial, U1Irrep)
+    Pspace = HO.hubbard_space(Trivial, U1Irrep)
     Vspace = Vect[FermionParity ⊠ U1Irrep]((0, 0) => 2, (1, 1 // 2) => 1, (1, -1 // 2) => 1)
     Espace = Vect[FermionParity ⊠ U1Irrep]((0, 0) => 8, (1, 1 // 2) => 4, (1, -1 // 2) => 4)
     truncs_env = collect(truncerror(; atol = 1.0e-12) & truncrank(χ) for χ in [8, 16])
diff --git a/test/timeevol/j1j2_finiteT.jl b/test/timeevol/j1j2_finiteT.jl
index 140fe6c19..4a590fe41 100644
--- a/test/timeevol/j1j2_finiteT.jl
+++ b/test/timeevol/j1j2_finiteT.jl
@@ -1,7 +1,6 @@
 using Test
 using LinearAlgebra
 using TensorKit
-import MPSKitModels: σˣ, σᶻ
 using PEPSKit
 
 # Benchmark energy from high-temperature expansion
diff --git a/test/timeevol/tf_ising_finiteT.jl b/test/timeevol/tf_ising_finiteT.jl
index e1137cc04..0845f3fba 100644
--- a/test/timeevol/tf_ising_finiteT.jl
+++ b/test/timeevol/tf_ising_finiteT.jl
@@ -1,7 +1,7 @@
 using Test
 using LinearAlgebra
 using TensorKit
-import MPSKitModels: σˣ, σᶻ
+import TensorKitTensors.SpinOperators as SO
 using PEPSKit
 
 # Benchmark data of [σx, σz] from HOTRG
@@ -21,8 +21,8 @@ end
 function measure_mag(pepo::InfinitePEPO, env::CTMRGEnv; purified::Bool = false)
     r, c = 1, 1
     lattice = physicalspace(pepo)
-    Mx = LocalOperator(lattice, ((r, c),) => σˣ(Float64, Trivial))
-    Mz = LocalOperator(lattice, ((r, c),) => σᶻ(Float64, Trivial))
+    Mx = LocalOperator(lattice, ((r, c),) => SO.σˣ(Float64, Trivial))
+    Mz = LocalOperator(lattice, ((r, c),) => SO.σᶻ(Float64, Trivial))
     if purified
         magx = expectation_value(pepo, Mx, pepo, env)
         magz = expectation_value(pepo, Mz, pepo, env)
diff --git a/test/types/localoperator.jl b/test/types/localoperator.jl
index dfca51c2b..eaf7b4049 100644
--- a/test/types/localoperator.jl
+++ b/test/types/localoperator.jl
@@ -2,7 +2,8 @@ using TensorKit
 using PEPSKit
 using PEPSKit: siterotl90, siterotr90, siterot180
 using MPSKit: add_physical_charge
-using MPSKitModels: a_number, nꜛnꜜ, contract_onesite
+import TensorKitTensors.BosonOperators as BO
+import TensorKitTensors.HubbardOperators as HO
 using Test
 
 vds = (ℂ^2, Rep[U₁](1 => 1, -1 => 1), Rep[SU₂](1 / 2 => 1))
@@ -54,8 +55,8 @@ end
 
     # bosonic case
     cutoff = 2
-    N = a_number(elt, U1Irrep; cutoff)
-    H_U = U / 2 * contract_onesite(N, N - id(domain(N)))
+    N = BO.b_num(elt, U1Irrep; cutoff)
+    H_U = U / 2 * N * (N - id(domain(N)))
     spaces = fill(space(H_U, 1), (lattice.Nrows, lattice.Ncols))
     H = LocalOperator(spaces, ((1, 1),) => H_U)
     tr_before = tr(last(only(H.terms)))
@@ -72,7 +73,7 @@ end
 
     # fermionic case
     symmetry = FermionParity ⊠ U1Irrep
-    H_U = U * nꜛnꜜ(elt, U1Irrep, Trivial)
+    H_U = U * HO.ud_num(elt, U1Irrep, Trivial)
     spaces = fill(space(H_U, 1), (lattice.Nrows, lattice.Ncols))
     H = LocalOperator(spaces, ((1, 1),) => H_U)
     tr_before = tr(last(only(H.terms)))
@@ -102,8 +103,10 @@ unitcells = [(1, 1), (2, 3), (3, 3), (4, 3)]
 end
 
 op_1x1 = LocalOperator([ℂ^2;;], ((1, 1), (1, 2)) => randn(ℂ^2, ℂ^2) ⊗ randn(ℂ^2, ℂ^2))
+# J1-J2 only has spin U(1) symmetry without sub-lattice rotation
+# See https://github.com/QuantumKitHub/MPSKitModels.jl/issues/57
 op_2x2 = add_physical_charge(
-    j1_j2_model(ComplexF64, U1Irrep, InfiniteSquare(2, 2)),
+    j1_j2_model(ComplexF64, U1Irrep, InfiniteSquare(2, 2); sublattice = false),
     [
         U1Irrep(-1 // 2) U1Irrep(1 // 2)
         U1Irrep(1 // 2) U1Irrep(-1 // 2)