From 02c06ce4e28dc5bbae1317e93a2567620ab5c34a Mon Sep 17 00:00:00 2001
From: Yue Zhengyuan <yuezy1997@icloud.com>
Date: Sun, 14 Jun 2026 16:21:40 +0800
Subject: [PATCH 1/9] Remove MPSKitModels dependency

- Replace all MPSKitModels operator building blocks with TensorKitTensors counterparts
  (SpinOperators, BosonOperators, HubbardOperators, TJOperators)
- Convert model constructors from MPSKitModels extensions to PEPSKit-native functions
- Remove AbstractLattice supertype; InfiniteSquare is now a plain struct
- Update tests, examples, docs, and notebooks to use TensorKitTensors imports
- Remove MPSKitModels from all Project.toml files

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
---
 Project.toml                               |  4 +-
 docs/Project.toml                          |  2 +-
 docs/make.jl                               |  6 +--
 docs/src/examples/bose_hubbard/index.md    |  3 +-
 docs/src/examples/bose_hubbard/main.ipynb  | 12 +----
 docs/src/examples/c4v_ctmrg/index.md       |  8 ++--
 docs/src/examples/c4v_ctmrg/main.ipynb     | 23 +--------
 docs/src/examples/heisenberg/index.md      |  3 +-
 docs/src/examples/heisenberg/main.ipynb    |  8 +---
 docs/src/examples/heisenberg_su/index.md   |  6 +--
 docs/src/examples/heisenberg_su/main.ipynb | 33 +------------
 docs/src/man/models.md                     | 55 ++++++++++++----------
 docs/src/man/precompilation.md             |  2 +-
 examples/Project.toml                      |  2 +-
 examples/bose_hubbard/main.jl              |  3 +-
 examples/c4v_ctmrg/main.jl                 |  8 ++--
 examples/heisenberg/main.jl                |  3 +-
 examples/heisenberg_su/main.jl             |  6 +--
 src/PEPSKit.jl                             |  6 ++-
 src/operators/lattices/squarelattice.jl    |  2 +-
 src/operators/models.jl                    | 51 ++++++++++----------
 test/Project.toml                          |  2 +-
 test/examples/heisenberg.jl                |  8 ++--
 test/timeevol/cluster_projectors.jl        |  4 +-
 test/timeevol/j1j2_finiteT.jl              |  2 +-
 test/timeevol/tf_ising_finiteT.jl          |  6 +--
 test/types/localoperator.jl                |  9 ++--
 27 files changed, 108 insertions(+), 169 deletions(-)

diff --git a/Project.toml b/Project.toml
index 3b36218a8..3e5cc7ba0 100644
--- a/Project.toml
+++ b/Project.toml
@@ -16,7 +16,6 @@ KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LoggingExtras = "e6f89c97-d47a-5376-807f-9c37f3926c36"
 MPSKit = "bb1c41ca-d63c-52ed-829e-0820dda26502"
-MPSKitModels = "ca635005-6f8c-4cd1-b51d-8491250ef2ab"
 MatrixAlgebraKit = "6c742aac-3347-4629-af66-fc926824e5e4"
 OhMyThreads = "67456a42-1dca-4109-a031-0a68de7e3ad5"
 OptimKit = "77e91f04-9b3b-57a6-a776-40b61faaebe0"
@@ -24,6 +23,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"
@@ -39,7 +39,6 @@ KrylovKit = "0.9.5, 0.10"
 LinearAlgebra = "1"
 LoggingExtras = "1"
 MPSKit = "0.13.9"
-MPSKitModels = "0.4"
 MatrixAlgebraKit = "0.6.5"
 OhMyThreads = "0.7, 0.8"
 OptimKit = "0.4"
@@ -47,6 +46,7 @@ Printf = "1"
 Random = "1"
 Statistics = "1"
 TensorKit = "0.16.5"
+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..75a16f5a4 100644
--- 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
index 636330ec3..249bec763 100644
--- 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
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..2ddcbe5b5 100644
--- a/docs/src/man/models.md
+++ b/docs/src/man/models.md
@@ -1,29 +1,40 @@
 # 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 particular, we use the submodules [`SpinOperators`](@extref TensorKitTensors SpinOperators),
+[`BosonOperators`](@extref TensorKitTensors BosonOperators),
+[`FermionOperators`](@extref TensorKitTensors FermionOperators),
+[`HubbardOperators`](@extref TensorKitTensors HubbardOperators), and
+[`TJOperators`](@extref TensorKitTensors TJOperators) to construct the local interaction tensors.
+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(
+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!(4 * SpinOperators.S_z_S_z(T, S), -J)
+    X = rmul!(SpinOperators.σˣ(T, S), g * -J)
     spaces = fill(domain(X)[1], (lattice.Nrows, lattice.Ncols))
     return LocalOperator(
         spaces,
@@ -35,7 +46,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`](@extref) 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,22 +54,16 @@ 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:
-
-### 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
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..58eedd2cc 100644
--- a/examples/Project.toml
+++ b/examples/Project.toml
@@ -3,7 +3,7 @@ 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"
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 66ae6e04e..7685f0f3e 100644
--- a/src/PEPSKit.jl
+++ b/src/PEPSKit.jl
@@ -33,7 +33,11 @@ using MPSKit: InfiniteEnvironments
 import MPSKit: tensorexpr, leading_boundary, loginit!, logiter!, logfinish!, logcancel!, physicalspace
 import MPSKit: infinite_temperature_density_matrix
 
-using MPSKitModels
+import TensorKitTensors.SpinOperators as SO
+import TensorKitTensors.FermionOperators as FO
+import TensorKitTensors.HubbardOperators as Hub
+import TensorKitTensors.TJOperators as tJ
+import TensorKitTensors.BosonOperators as BO
 using FiniteDifferences
 using OhMyThreads: tmap, tmap!
 using DocStringExtensions
diff --git a/src/operators/lattices/squarelattice.jl b/src/operators/lattices/squarelattice.jl
index 102546594..b42457f81 100644
--- a/src/operators/lattices/squarelattice.jl
+++ b/src/operators/lattices/squarelattice.jl
@@ -13,7 +13,7 @@ $(TYPEDFIELDS)
 
 By default, an infinite square with a (1, 1)-unitcell is constructed.
 """
-struct InfiniteSquare <: AbstractLattice{2}
+struct InfiniteSquare
     Nrows::Int
     Ncols::Int
     function InfiniteSquare(Nrows::Integer = 1, Ncols::Integer = 1)
diff --git a/src/operators/models.jl b/src/operators/models.jl
index 6e1d927b9..eb5d8dc42 100644
--- a/src/operators/models.jl
+++ b/src/operators/models.jl
@@ -22,15 +22,15 @@ function nearest_neighbour_hamiltonian(
 end
 
 #
-## MPSKitModels.jl re-exports
+## Model definitions
 #
 
-function MPSKitModels.transverse_field_ising(
+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!(4 * SO.S_z_S_z(T, S), -J)
+    X = rmul!(SO.σˣ(T, S), g * -J)
     spaces = fill(domain(X)[1], (lattice.Nrows, lattice.Ncols))
     return LocalOperator(
         spaces,
@@ -39,31 +39,31 @@ function MPSKitModels.transverse_field_ising(
     )
 end
 
-function MPSKitModels.heisenberg_XYZ(lattice::InfiniteSquare; kwargs...)
+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(
+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 +71,34 @@ function MPSKitModels.heisenberg_XXZ(
     )
 end
 
-function MPSKitModels.hubbard_model(
+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 = Hub.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)
+        Hub.e_plus_e_min(T, particle_symmetry, spin_symmetry) +
+        Hub.e_min_e_plus(T, particle_symmetry, spin_symmetry)
+    interaction_term = Hub.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(
+function bose_hubbard_model(
         elt::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)))
+        BO.b_plus_b_min(elt, symmetry; cutoff = cutoff) +
+        BO.b_min_b_plus(elt, symmetry; cutoff = cutoff)
+    N = BO.b_num(elt, symmetry; cutoff = cutoff)
+    interaction_term = N * (N - id(domain(N)))
 
     spaces = fill(space(N, 1), (lattice.Nrows, lattice.Ncols))
 
@@ -120,7 +121,7 @@ function MPSKitModels.bose_hubbard_model(
     return H
 end
 
-function MPSKitModels.tj_model(
+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,
@@ -142,7 +143,7 @@ function MPSKitModels.tj_model(
 end
 
 #
-## Models not yet defined in MPSKitModels
+## Additional models
 #
 
 """
@@ -168,9 +169,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
diff --git a/test/Project.toml b/test/Project.toml
index 2742a2812..40fc24950 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -7,7 +7,7 @@ ChainRulesTestUtils = "cdddcdb0-9152-4a09-a978-84456f9df70a"
 KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MPSKit = "bb1c41ca-d63c-52ed-829e-0820dda26502"
-MPSKitModels = "ca635005-6f8c-4cd1-b51d-8491250ef2ab"
+
 MatrixAlgebraKit = "6c742aac-3347-4629-af66-fc926824e5e4"
 OptimKit = "77e91f04-9b3b-57a6-a776-40b61faaebe0"
 PEPSKit = "52969e89-939e-4361-9b68-9bc7cde4bdeb"
diff --git a/test/examples/heisenberg.jl b/test/examples/heisenberg.jl
index 6cb797b2b..8c018893b 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..7511ac073 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 Hub
 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 = Hub.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 51009cbe1..6975998c2 100644
--- a/test/timeevol/j1j2_finiteT.jl
+++ b/test/timeevol/j1j2_finiteT.jl
@@ -1,7 +1,7 @@
 using Test
 using LinearAlgebra
 using TensorKit
-import MPSKitModels: σˣ, σᶻ
+import TensorKitTensors.SpinOperators as SO
 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 62cc7edc7..6ac805ded 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..db9a5f98e 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 Hub
 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 * Hub.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)))

From 8ca4b37ba6efdfe3fe0dbf4afb07ae8b198cff06 Mon Sep 17 00:00:00 2001
From: Yue Zhengyuan <yuezy1997@icloud.com>
Date: Sun, 14 Jun 2026 16:39:30 +0800
Subject: [PATCH 2/9] Rewrite pwave_superconductor

---
 src/operators/models.jl | 27 +++++++++++++--------------
 1 file changed, 13 insertions(+), 14 deletions(-)

diff --git a/src/operators/models.jl b/src/operators/models.jl
index eb5d8dc42..d26d75649 100644
--- a/src/operators/models.jl
+++ b/src/operators/models.jl
@@ -186,15 +186,17 @@ end
 """
     pwave_superconductor([T=ComplexF64,] lattice::InfiniteSquare; t=1, μ=2, Δ=1)
 
-Square lattice ``p``-wave superconductor model, defined by the Hamiltonian
+Square lattice ``(p + ip)``-wave superconductor model, defined by the Hamiltonian
 
 ```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...)
@@ -203,22 +205,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]
+    hx = 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))

From 6deceb2cb447a88b8714d1108df0768ae6d350da Mon Sep 17 00:00:00 2001
From: Yue Zhengyuan <yuezy1997@icloud.com>
Date: Sun, 14 Jun 2026 17:01:29 +0800
Subject: [PATCH 3/9] Fix hopping terms

---
 src/operators/models.jl | 26 ++++++++------------------
 1 file changed, 8 insertions(+), 18 deletions(-)

diff --git a/src/operators/models.jl b/src/operators/models.jl
index d26d75649..259a8e280 100644
--- a/src/operators/models.jl
+++ b/src/operators/models.jl
@@ -81,9 +81,7 @@ function hubbard_model(
     N = Hub.e_num(T, particle_symmetry, spin_symmetry)
     pspace = space(N, 1)
     unit = TensorKit.id(pspace)
-    hopping =
-        Hub.e_plus_e_min(T, particle_symmetry, spin_symmetry) +
-        Hub.e_min_e_plus(T, particle_symmetry, spin_symmetry)
+    hopping = Hub.e_hopping(T, particle_symmetry, spin_symmetry)
     interaction_term = Hub.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)
@@ -91,13 +89,11 @@ function hubbard_model(
 end
 
 function bose_hubbard_model(
-        elt::Type{<:Number}, symmetry::Type{<:Sector}, lattice::InfiniteSquare;
+        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 =
-        BO.b_plus_b_min(elt, symmetry; cutoff = cutoff) +
-        BO.b_min_b_plus(elt, symmetry; cutoff = cutoff)
-    N = BO.b_num(elt, symmetry; cutoff = cutoff)
+    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))
@@ -126,19 +122,13 @@ function tj_model(
         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_number(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
 
@@ -147,7 +137,7 @@ end
 #
 
 """
-    j1_j2_model([elt::Type{T}, symm::Type{S},] lattice::InfiniteSquare;
+    j1_j2_model([T::Type{T}, symm::Type{S},] 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

From d95806e9d948b0b3c28eb4bf21b0ce81a73fb62f Mon Sep 17 00:00:00 2001
From: Yue Zhengyuan <yuezy1997@icloud.com>
Date: Sun, 14 Jun 2026 17:01:35 +0800
Subject: [PATCH 4/9] Port over some docstrings

---
 src/operators/lattices/squarelattice.jl | 15 ++++
 src/operators/models.jl                 | 91 ++++++++++++++++++++++++-
 2 files changed, 105 insertions(+), 1 deletion(-)

diff --git a/src/operators/lattices/squarelattice.jl b/src/operators/lattices/squarelattice.jl
index b42457f81..1b1242d32 100644
--- a/src/operators/lattices/squarelattice.jl
+++ b/src/operators/lattices/squarelattice.jl
@@ -24,10 +24,20 @@ end
 
 Base.size(lattice::InfiniteSquare) = (lattice.Nrows, lattice.Ncols)
 
+"""
+    vertices(lattice::InfiniteSquare)
+
+Return an iterator over all vertices in the unit cell.
+"""
 function vertices(lattice::InfiniteSquare)
     return CartesianIndices((1:(lattice.Nrows), 1:(lattice.Ncols)))
 end
 
+"""
+    nearest_neighbours(lattice::InfiniteSquare)
+
+Return an iterator over all pairs of nearest neighbours.
+"""
 function nearest_neighbours(lattice::InfiniteSquare)
     neighbors = Vector{CartesianIndex{2}}[]
     for idx in vertices(lattice)
@@ -37,6 +47,11 @@ function nearest_neighbours(lattice::InfiniteSquare)
     return neighbors
 end
 
+"""
+    next_nearest_neighbours(lattice::InfiniteSquare)
+
+Return an iterator over all pairs of next-nearest neighbours.
+"""
 function next_nearest_neighbours(lattice::InfiniteSquare)
     neighbors = Vector{CartesianIndex{2}}[]
     for idx in vertices(lattice)
diff --git a/src/operators/models.jl b/src/operators/models.jl
index 259a8e280..cb5cac194 100644
--- a/src/operators/models.jl
+++ b/src/operators/models.jl
@@ -25,6 +25,19 @@ end
 ## Model definitions
 #
 
+"""
+    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)
+on an infinite square lattice, as defined by
+```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, the model is defined with `Trivial` symmetry and with `ComplexF64` entries.
+"""
 function transverse_field_ising(
         T::Type{<:Number}, S::Union{Type{Trivial}, Type{Z2Irrep}}, lattice::InfiniteSquare;
         J = 1.0, g = 1.0,
@@ -39,6 +52,21 @@ function transverse_field_ising(
     )
 end
 
+"""
+    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 on an infinite square lattice, as defined by
+```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, the model uses the antiferromagnetic convention ``(J_x, J_y, J_z) = (-1, 1, -1)``
+which is suitable for a single-site unit cell ground state after a sublattice rotation.
+It is defined with `Trivial` symmetry and with `ComplexF64` entries.
+
+See also [`heisenberg_XXZ`](@ref).
+"""
 function heisenberg_XYZ(lattice::InfiniteSquare; kwargs...)
     return heisenberg_XYZ(ComplexF64, Trivial, lattice; kwargs...)
 end
@@ -56,6 +84,19 @@ function heisenberg_XYZ(
     )
 end
 
+"""
+    heisenberg_XXZ([T::Type{<:Number}], [S::Type{<:Sector}], lattice::InfiniteSquare;
+                   J=1.0, Delta=1.0, spin=1)
+
+`LocalOperator` for the XXZ Heisenberg model on an infinite square lattice, as defined by
+```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, the model is defined with `Trivial` symmetry and with `ComplexF64` entries.
+
+See also [`heisenberg_XYZ`](@ref).
+"""
 function heisenberg_XXZ(
         T::Type{<:Number}, S::Type{<:Sector}, lattice::InfiniteSquare;
         J = 1.0, Delta = 1.0, spin = 1
@@ -71,6 +112,21 @@ function heisenberg_XXZ(
     )
 end
 
+"""
+    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 Hubbard model on an infinite square lattice, as defined by
+```math
+H = -t \\sum_{\\langle i,j \\rangle} \\sum_{\\sigma} \\left( e_{i,\\sigma}^\\dagger e_{j,\\sigma} + e_{j,\\sigma}^\\dagger e_{i,\\sigma} \\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, the model is defined without any symmetries and with `ComplexF64` entries.
+"""
 function hubbard_model(
         T::Type{<:Number}, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector},
         lattice::InfiniteSquare;
@@ -88,6 +144,22 @@ function hubbard_model(
     return nearest_neighbour_hamiltonian(fill(pspace, size(lattice)), h)
 end
 
+"""
+    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 on an infinite square lattice, as defined by
+```math
+H = -t \\sum_{\\langle i,j \\rangle} \\left( a_i^\\dagger a_j + a_j^\\dagger a_i \\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, the model is defined with `Trivial` symmetry and `ComplexF64` entries.
+"""
 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,
@@ -117,6 +189,23 @@ function bose_hubbard_model(
     return H
 end
 
+"""
+    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 on an infinite square lattice, as defined by
+```math
+H = -t \\sum_{\\langle i,j \\rangle, \\sigma}
+    (\\tilde{e}^\\dagger_{i,\\sigma} \\tilde{e}_{j,\\sigma} + 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`` restricted to
+the no-double-occupancy subspace.
+
+By default, the model is defined without any symmetries and with `ComplexF64` entries.
+"""
 function tj_model(
         T::Type{<:Number}, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector},
         lattice::InfiniteSquare;
@@ -207,7 +296,7 @@ function pwave_superconductor(
     # two-site (x-direction)
     hx = hopp + (pair_x + pair_x')
     # two-site (y-direction)
-    hx = hopp + (pair_y + pair_y')
+    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))

From f4418136530b780f0ae0d7c87d381ea1e087c7f1 Mon Sep 17 00:00:00 2001
From: Yue Zhengyuan <yuezy1997@icloud.com>
Date: Sun, 14 Jun 2026 17:13:37 +0800
Subject: [PATCH 5/9] Update Project.toml

---
 examples/Project.toml         | 2 +-
 test/Project.toml             | 2 +-
 test/timeevol/j1j2_finiteT.jl | 1 -
 3 files changed, 2 insertions(+), 3 deletions(-)

diff --git a/examples/Project.toml b/examples/Project.toml
index 58eedd2cc..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"
-
 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/test/Project.toml b/test/Project.toml
index 40fc24950..56cfe7eaf 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -7,7 +7,6 @@ ChainRulesTestUtils = "cdddcdb0-9152-4a09-a978-84456f9df70a"
 KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MPSKit = "bb1c41ca-d63c-52ed-829e-0820dda26502"
-
 MatrixAlgebraKit = "6c742aac-3347-4629-af66-fc926824e5e4"
 OptimKit = "77e91f04-9b3b-57a6-a776-40b61faaebe0"
 PEPSKit = "52969e89-939e-4361-9b68-9bc7cde4bdeb"
@@ -15,6 +14,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/timeevol/j1j2_finiteT.jl b/test/timeevol/j1j2_finiteT.jl
index 6975998c2..7cf888fda 100644
--- a/test/timeevol/j1j2_finiteT.jl
+++ b/test/timeevol/j1j2_finiteT.jl
@@ -1,7 +1,6 @@
 using Test
 using LinearAlgebra
 using TensorKit
-import TensorKitTensors.SpinOperators as SO
 using PEPSKit
 
 # Benchmark energy from high-temperature expansion

From 4d946144126713ce44442438d4530c7324286e54 Mon Sep 17 00:00:00 2001
From: Yue Zhengyuan <yuezy1997@icloud.com>
Date: Sun, 14 Jun 2026 17:17:44 +0800
Subject: [PATCH 6/9] Remove broken doc links

---
 docs/src/man/models.md | 7 +------
 1 file changed, 1 insertion(+), 6 deletions(-)

diff --git a/docs/src/man/models.md b/docs/src/man/models.md
index 2ddcbe5b5..d75fe1482 100644
--- a/docs/src/man/models.md
+++ b/docs/src/man/models.md
@@ -9,11 +9,6 @@ already implemented models.
 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 particular, we use the submodules [`SpinOperators`](@extref TensorKitTensors SpinOperators),
-[`BosonOperators`](@extref TensorKitTensors BosonOperators),
-[`FermionOperators`](@extref TensorKitTensors FermionOperators),
-[`HubbardOperators`](@extref TensorKitTensors HubbardOperators), and
-[`TJOperators`](@extref TensorKitTensors TJOperators) to construct the local interaction tensors.
 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),
@@ -46,7 +41,7 @@ end
 
 This provides a good recipe for defining a model:
 
-1. Define the locally-acting tensors as `TensorMap`s using [`TensorKitTensors`](@extref) operators.
+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.
 

From 2143b25834c22c5a46cb864b289dc2baebe9b3f6 Mon Sep 17 00:00:00 2001
From: Yue Zhengyuan <yuezy1997@icloud.com>
Date: Sun, 14 Jun 2026 17:36:43 +0800
Subject: [PATCH 7/9] Fix add_physical_charge test

---
 test/types/localoperator.jl | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/test/types/localoperator.jl b/test/types/localoperator.jl
index db9a5f98e..4f0dcaab7 100644
--- a/test/types/localoperator.jl
+++ b/test/types/localoperator.jl
@@ -103,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)

From 866d20737a298e3ebcd7433849611e3828960acf Mon Sep 17 00:00:00 2001
From: Yue Zhengyuan <yuezy1997@icloud.com>
Date: Sun, 14 Jun 2026 18:44:51 +0800
Subject: [PATCH 8/9] Restore default T and sym methods for Hamiltonians

---
 src/operators/models.jl | 17 ++++++++++++++++-
 1 file changed, 16 insertions(+), 1 deletion(-)

diff --git a/src/operators/models.jl b/src/operators/models.jl
index cb5cac194..0f4b1f6ee 100644
--- a/src/operators/models.jl
+++ b/src/operators/models.jl
@@ -38,6 +38,9 @@ where ``\\sigma^i`` are the spin-1/2 Pauli operators.
 
 By default, the model is defined with `Trivial` symmetry and with `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,
@@ -97,6 +100,9 @@ By default, the model is defined with `Trivial` symmetry and with `ComplexF64` e
 
 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
@@ -127,6 +133,9 @@ fermionic number operator.
 
 By default, the model 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;
@@ -160,6 +169,9 @@ If `symmetry` is `U1Irrep`, a fixed (half-integer) particle number density `n` c
 
 By default, the model 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,
@@ -206,13 +218,16 @@ the no-double-occupancy subspace.
 
 By default, the model 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 = tJ.e_hopping(T, particle_symmetry, spin_symmetry; slave_fermion)
-    num = tJ.e_number(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)

From ccf00ffb888a83441eaa5138cf4423d34953e895 Mon Sep 17 00:00:00 2001
From: Yue Zhengyuan <yuezy1997@icloud.com>
Date: Sun, 14 Jun 2026 19:00:24 +0800
Subject: [PATCH 9/9] Change p-wave convention

---
 src/operators/models.jl | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/operators/models.jl b/src/operators/models.jl
index 0f4b1f6ee..d599025e3 100644
--- a/src/operators/models.jl
+++ b/src/operators/models.jl
@@ -280,7 +280,7 @@ end
 """
     pwave_superconductor([T=ComplexF64,] lattice::InfiniteSquare; t=1, μ=2, Δ=1)
 
-Square lattice ``(p + ip)``-wave superconductor model, defined by the Hamiltonian
+Square lattice ``(p - ip)``-wave superconductor model, defined by the Hamiltonian
 
 ```math
     H = -\\sum_{\\langle i,j \\rangle} \\left( t c_i^\\dagger c_j +
@@ -289,8 +289,8 @@ Square lattice ``(p + ip)``-wave superconductor model, defined by the Hamiltonia
 
 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.
+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...)
@@ -304,7 +304,7 @@ function pwave_superconductor(
     # 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
+    pair_x, pair_y = Δ * pair, -im * Δ * pair
 
     # on-site
     h0 = -μ * FO.f_num(T, Trivial)