Yue-Zhengyuan · Aug 2, 2026 · Jun 14, 2026 · Jun 14, 2026 · Jun 14, 2026 · Jun 14, 2026
diff --git 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.3"
 TensorOperations = "5"
 TupleTools = "1.6.0"
 VectorInterface = "0.4, 0.5, 0.6"

diff --git 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
@@ -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
@@ -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
@@ -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
@@ -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
+using TensorKitTensors.SpinOperators
 
 # 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
+        S_x_S_x(T, S; spin = spin) * Jx +
+        S_y_S_y(T, S; spin = spin) * Jy +
+        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
@@ -62,7 +62,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "using MPSKitModels: S_xx, S_yy, S_zz\n",
+    "using TensorKitTensors.SpinOperators\n",
     "\n",
     "# Heisenberg model assuming C4v symmetric PEPS and environment, which only evaluates necessary term\n",
     "function heisenberg_XYZ_c4v(lattice::InfiniteSquare; kwargs...)\n",
@@ -74,9 +74,9 @@
     "    )\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",
+    "        S_x_S_x(T, S; spin = spin) * Jx +\n",
+    "        S_y_S_y(T, S; spin = spin) * Jy +\n",
+    "        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",
@@ -293,4 +293,4 @@
  },
  "nbformat": 4,
  "nbformat_minor": 3
-}
+}
diff --git 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
@@ -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
@@ -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
+using TensorKitTensors.SpinOperators
 Random.seed!(0);
 ````
 

diff --git a/docs/src/examples/heisenberg_su/main.ipynb b/docs/src/examples/heisenberg_su/main.ipynb
@@ -38,7 +38,7 @@
     "using Random\n",
     "import Statistics: mean\n",
     "using TensorKit, PEPSKit\n",
-    "import MPSKitModels: S_x, S_y, S_z, S_exchange\n",
+    "using TensorKitTensors.SpinOperators\n",
     "Random.seed!(0);"
    ]
   },
@@ -249,4 +249,4 @@
  },
  "nbformat": 4,
  "nbformat_minor": 3
-}
+}
diff --git a/docs/src/man/models.md b/docs/src/man/models.md
@@ -1,28 +1,31 @@
 # 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(
+using TensorKitTensors.SpinOperators
+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)
+    ZZ = rmul!(S_z_S_z(T, S), -4 * J)
     X = rmul!(σˣ(T, S), g * -J)
     spaces = fill(domain(X)[1], (lattice.Nrows, lattice.Ncols))
     return LocalOperator(
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
+using TensorKitTensors.SpinOperators
 
 ## 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
+        S_x_S_x(T, S; spin = spin) * Jx +
+        S_y_S_y(T, S; spin = spin) * Jy +
+        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
@@ -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
@@ -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
+using TensorKitTensors.SpinOperators
 Random.seed!(0);
 
 md"""

diff --git 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