QuantumKitHub · Yue-Zhengyuan · 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.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
@@ -4,7 +4,7 @@ DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244"
 DocumenterInterLinks = "d12716ef-a0f6-4df4-a9f1-a5a34e75c656"
 KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77"
 MPSKit = "bb1c41ca-d63c-52ed-829e-0820dda26502"
-MPSKitModels = "ca635005-6f8c-4cd1-b51d-8491250ef2ab"
+TensorKitTensors = "41b62e7d-e9d1-4e23-942c-79a97adf954b"
 PEPSKit = "52969e89-939e-4361-9b68-9bc7cde4bdeb"
 TensorKit = "07d1fe3e-3e46-537d-9eac-e9e13d0d4cec"
 

diff --git a/docs/make.jl b/docs/make.jl
@@ -11,7 +11,6 @@ using Documenter
 using DocumenterCitations
 using DocumenterInterLinks
 using PEPSKit
-using MPSKitModels: MPSKitModels # used for docstrings
 
 # bibliography
 bibpath = joinpath(@__DIR__, "src", "assets", "pepskit.bib")
@@ -23,7 +22,7 @@ links = InterLinks(
     "TensorKit" => "https://quantumkithub.github.io/TensorKit.jl/stable/",
     "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/",
@@ -54,7 +53,7 @@ examples_partition_functions = joinpath.(
 examples_boundary_mps = joinpath.(["boundary_mps"], Ref("index.md"))
 
 makedocs(;
-    modules = [PEPSKit, MPSKitModels],
+    modules = [PEPSKit],
     sitename = "PEPSKit.jl",
     format = Documenter.HTML(;
         prettyurls = true, mathengine, assets = ["assets/custom.css"], size_threshold = 1024000
@@ -75,7 +74,6 @@ makedocs(;
         "References" => "references.md",
     ],
     checkdocs = :none,
-    # checkdocs_ignored_modules=[MPSKitModels], # doesn't seem to work...
     plugins = [bib, links],
 )
 

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
+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
@@ -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
@@ -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
+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
@@ -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
@@ -1,29 +1,36 @@
 # 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 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).
+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 lattice types such as `InfiniteSquare`.
+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,30 +42,24 @@ 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.
 
 For more model implementations, check the [PEPSKit repository](https://github.com/QuantumKitHub/PEPSKit.jl/blob/master/src/operators/models.jl).
 
 ## 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:
-
-### MPSKitModels.jl models
-
-```@docs
-MPSKitModels.transverse_field_ising
-MPSKitModels.heisenberg_XYZ
-MPSKitModels.heisenberg_XXZ
-MPSKitModels.hubbard_model
-MPSKitModels.bose_hubbard_model
-MPSKitModels.tj_model
-```
-
-### PEPSKit.jl models
+PEPSKit provides a number of pre-defined models. The following model constructors are available
+and can be used directly with an [`InfiniteSquare`](@ref) lattice:
 
 ```@docs
+transverse_field_ising
+heisenberg_XYZ
+heisenberg_XXZ
+hubbard_model
+bose_hubbard_model
+tj_model
 j1_j2_model
 pwave_superconductor
 ```
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
+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
@@ -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)