Add t-J model

.gitignore

@@ -1,2 +1,4 @@
 docs/build
 Manifest.toml
+.vscode/
+.DS_Store

docs/make.jl

@@ -9,10 +9,8 @@ makedocs(;
                                 prettyurls=get(ENV, "CI", nothing) == "true",
                                 mathengine=MathJax()),
          pages=["Home" => "index.md",
-                "Manual" => ["man/operators.md",
-                             "man/mpoham.md",
-                             "man/lattices.md",
+                "Manual" => ["man/operators.md", "man/mpoham.md", "man/lattices.md",
                              "man/models.md"],
-                "Index" => "package_index.md"])
+                "Index" => "package_index.md"],)
 
 deploydocs(; repo="github.com/QuantumKitHub/MPSKitModels.jl.git")

src/MPSKitModels.jl

@@ -35,13 +35,15 @@ export c⁺, c⁻, c⁺⁺, c⁻⁻, c⁺⁻, c⁻⁺
 export e_plus, e_min, e_plusplus, e_minmin, e_plusmin, e_minplus
 export e_number, e_number_up, e_number_down, e_number_updown
 export e⁺⁺, e⁻⁻, e⁺⁻, e⁻⁺
+export tJ, TJOperators
 
 export transverse_field_ising
 export kitaev_model
 export quantum_potts
 export heisenberg_XXX, heisenberg_XXZ, heisenberg_XYZ
 export bilinear_biquadratic_model
 export hubbard_model, bose_hubbard_model
+export tj_model
 export quantum_chemistry_hamiltonian
 
 export classical_ising
@@ -65,6 +67,11 @@ include("operators/spinoperators.jl")
 include("operators/fermionoperators.jl")
 include("operators/hubbardoperators.jl")
 using .HubbardOperators
+# TJOperators share operator names with HubbardOperators
+# and is only imported to avoid name conflicts
+include("operators/tjoperators.jl")
+import .TJOperators
+const tJ = TJOperators
 include("operators/bosonoperators.jl")
 
 include("models/hamiltonians.jl")

src/lattices/squarelattice.jl

@@ -149,28 +149,36 @@
 function nearest_neighbours(lattice::FiniteStrip)
     rows = lattice.L
     cols = lattice.N ÷ lattice.L
-    horizontal = (LatticePoint((i, j), lattice) => LatticePoint((i, j + 1), lattice)
-                  for i in 1:rows, j in 1:(cols - 1))
-    vertical = (LatticePoint((i, j), lattice) => LatticePoint((i + 1, j), lattice)
-                for i in 1:(rows - 1), j in 1:cols)
+    horizontal = (LatticePoint((i, j), lattice) => LatticePoint((i, j + 1), lattice) for i in
+                                                                                         1:rows,
+                                                                                         j in
+                                                                                         1:(cols - 1))
+    vertical = (LatticePoint((i, j), lattice) => LatticePoint((i + 1, j), lattice) for
+                i in 1:(rows - 1), j in 1:cols)
     return Iterators.flatten((horizontal, vertical))
 end
 function nearest_neighbours(lattice::FiniteCylinder)
     rows = lattice.L
     cols = lattice.N ÷ lattice.L
-    horizontal = (LatticePoint((i, j), lattice) => LatticePoint((i, j + 1), lattice)
-                  for i in 1:rows, j in 1:(cols - 1))
-    vertical = (LatticePoint((i, j), lattice) => LatticePoint((i + 1, j), lattice)
-                for i in 1:rows, j in 1:cols)
+    horizontal = (LatticePoint((i, j), lattice) => LatticePoint((i, j + 1), lattice) for i in
+                                                                                         1:rows,
+                                                                                         j in
+                                                                                         1:(cols - 1))
+    vertical = (LatticePoint((i, j), lattice) => LatticePoint((i + 1, j), lattice) for i in
+                                                                                       1:rows,
+                                                                                       j in
+                                                                                       1:cols)
     return Iterators.flatten((horizontal, vertical))
 end
 function nearest_neighbours(lattice::FiniteHelix)
     rows = lattice.L
     cols = lattice.N ÷ lattice.L
-    horizontal = (LatticePoint((i, j), lattice) => LatticePoint((i, j + 1), lattice)
-                  for i in 1:rows, j in 1:(cols - 1))
-    vertical = (LatticePoint((i, j), lattice) => LatticePoint((i + 1, j), lattice)
-                for i in 1:rows, j in 1:cols if (i != rows && j != cols))
+    horizontal = (LatticePoint((i, j), lattice) => LatticePoint((i, j + 1), lattice) for i in
+                                                                                         1:rows,
+                                                                                         j in
+                                                                                         1:(cols - 1))
+    vertical = (LatticePoint((i, j), lattice) => LatticePoint((i + 1, j), lattice) for
+                i in 1:rows, j in 1:cols if (i != rows && j != cols))
     return Iterators.flatten((horizontal, vertical))
 end
 function nearest_neighbours(lattice::Union{InfiniteStrip,InfiniteCylinder,InfiniteHelix})

src/models/hamiltonians.jl

@@ -23,23 +23,21 @@
     return transverse_field_ising(ComplexF64, Trivial, lattice; kwargs...)
 end
 function transverse_field_ising(S::Type{<:Sector},
-                                lattice::AbstractLattice=InfiniteChain(1);
-                                kwargs...)
+                                lattice::AbstractLattice=InfiniteChain(1); kwargs...)
     return transverse_field_ising(ComplexF64, S, lattice; kwargs...)
 end
 function transverse_field_ising(T::Type{<:Number}, lattice::AbstractLattice; kwargs...)
     return transverse_field_ising(T, Trivial, lattice; kwargs...)
 end
-function transverse_field_ising(T::Type{<:Number},
-                                S::Type{<:Sector},
-                                lattice::AbstractLattice;
-                                kwargs...)
+function transverse_field_ising(T::Type{<:Number}, S::Type{<:Sector},
+                                lattice::AbstractLattice; kwargs...)
     throw(ArgumentError("`symmetry` must be either `Trivial`, `Z2Irrep` or `FermionParity`"))
 end
 function transverse_field_ising(T::Type{<:Number}=ComplexF64,
                                 S::Union{Type{Trivial},Type{Z2Irrep}}=Trivial,
                                 lattice::AbstractLattice=InfiniteChain(1);
-                                J=1.0, g=1.0)
+                                J=1.0,
+                                g=1.0,)
     ZZ = rmul!(σᶻᶻ(T, S), -J)
     X = rmul!(σˣ(T, S), g * -J)
     return @mpoham begin
@@ -51,8 +49,7 @@
     end
 end
 function transverse_field_ising(T::Type{<:Number}, ::Type{fℤ₂},
-                                lattice::AbstractLattice=InfiniteChain(1);
-                                J=1.0, g=1.0)
+                                lattice::AbstractLattice=InfiniteChain(1); J=1.0, g=1.0)
     twosite = axpby!(-J, c_plusmin(T) + c_minplus(T), J, c_plusplus(T) + c_minmin(T))
     onesite = axpby!(2g * J, c_number(T), -g * J, id(Matrix{T}, space(twosite, 1)))
 
@@ -86,7 +83,9 @@
 end
 function kitaev_model(elt::Type{<:Number}=ComplexF64,
                       lattice::AbstractLattice=InfiniteChain(1);
-                      t=1.0, mu=1.0, Delta=1.0)
+                      t=1.0,
+                      mu=1.0,
+                      Delta=1.0,)
     TB = rmul!(c_plusmin(elt) + c_minplus(elt), -t / 2)     # tight-binding term
     SC = rmul!(c_plusplus(elt) + c_minmin(elt), Delta / 2)  # superconducting term
     CP = rmul!(c_number(elt), -mu)                          # chemical potential term
@@ -132,7 +131,8 @@
 function heisenberg_XXX(T::Type{<:Number}=ComplexF64,
                         symmetry::Type{<:Sector}=Trivial,
                         lattice::AbstractLattice=InfiniteChain(1);
-                        J::Real=1.0, spin::Real=1)
+                        J::Real=1.0,
+                        spin::Real=1,)
     term = rmul!(S_exchange(T, symmetry; spin=spin), J)
     return @mpoham sum(nearest_neighbours(lattice)) do (i, j)
         return term{i,j}
@@ -164,7 +164,9 @@
 function heisenberg_XXZ(elt::Type{<:Number}=ComplexF64,
                         symmetry::Type{<:Sector}=Trivial,
                         lattice::AbstractLattice=InfiniteChain(1);
-                        J=1.0, Delta=1.0, spin=1)
+                        J=1.0,
+                        Delta=1.0,
+                        spin=1,)
     term = rmul!(S_xx(elt, symmetry; spin=spin), J) +
            rmul!(S_yy(elt, symmetry; spin=spin), J) +
            rmul!(S_zz(elt, symmetry; spin=spin), Delta * J)
@@ -190,7 +192,10 @@
 end
 function heisenberg_XYZ(T::Type{<:Number}=ComplexF64,
                         lattice::AbstractLattice=InfiniteChain(1);
-                        Jx=1.0, Jy=1.0, Jz=1.0, spin=1)
+                        Jx=1.0,
+                        Jy=1.0,
+                        Jz=1.0,
+                        spin=1,)
     term = rmul!(S_xx(T, Trivial; spin=spin), Jx) +
            rmul!(S_yy(T, Trivial; spin=spin), Jy) +
            rmul!(S_zz(T, Trivial; spin=spin), Jz)
@@ -226,7 +231,9 @@
 function bilinear_biquadratic_model(elt::Type{<:Number}=ComplexF64,
                                     symmetry::Type{<:Sector}=Trivial,
                                     lattice::AbstractLattice=InfiniteChain(1);
-                                    spin=1, J=1.0, θ=0.0)
+                                    spin=1,
+                                    J=1.0,
+                                    θ=0.0,)
     return @mpoham sum(nearest_neighbours(lattice)) do (i, j)
         return J * cos(θ) * S_exchange(elt, symmetry; spin=spin){i,j} +
                J * sin(θ) * (S_exchange(elt, symmetry; spin=spin)^2){i,j}
@@ -252,18 +259,19 @@
 function quantum_potts(lattice::AbstractLattice; kwargs...)
     return quantum_potts(ComplexF64, Trivial, lattice; kwargs...)
 end
-function quantum_potts(symmetry::Type{<:Sector},
-                       lattice::AbstractLattice=InfiniteChain(1); kwargs...)
+function quantum_potts(symmetry::Type{<:Sector}, lattice::AbstractLattice=InfiniteChain(1);
+                       kwargs...)
     return quantum_potts(ComplexF64, symmetry, lattice; kwargs...)
 end
-function quantum_potts(elt::Type{<:Number}, lattice::AbstractLattice;
-                       kwargs...)
+function quantum_potts(elt::Type{<:Number}, lattice::AbstractLattice; kwargs...)
     return quantum_potts(elt, Trivial, lattice; kwargs...)
 end
 function quantum_potts(elt::Type{<:Number}=ComplexF64,
                        symmetry::Type{<:Sector}=Trivial,
                        lattice::AbstractLattice=InfiniteChain(1);
-                       q=3, J=1.0, g=1.0)
+                       q=3,
+                       J=1.0,
+                       g=1.0,)
     return @mpoham sum(sum(nearest_neighbours(lattice)) do (i, j)
                            return -J * (potts_ZZ(elt, symmetry; q)^k){i,j}
                        end - sum(vertices(lattice)) do i
@@ -304,18 +312,20 @@
                        particle_symmetry::Type{<:Sector}=Trivial,
                        spin_symmetry::Type{<:Sector}=Trivial,
                        lattice::AbstractLattice=InfiniteChain(1);
-                       t=1.0, U=1.0, mu=0.0, n::Integer=0)
+                       t=1.0,
+                       U=1.0,
+                       mu=0.0,
+                       n::Integer=0,)
     hopping = e⁺e⁻(T, particle_symmetry, spin_symmetry) +
               e⁻e⁺(T, particle_symmetry, spin_symmetry)
     interaction_term = nꜛnꜜ(T, particle_symmetry, spin_symmetry)
     N = e_number(T, particle_symmetry, spin_symmetry)
     return @mpoham begin
         sum(nearest_neighbours(lattice)) do (i, j)
             return -t * hopping{i,j}
-        end +
-        sum(vertices(lattice)) do i
-            return U * interaction_term{i} - mu * N{i}
-        end
+        end + sum(vertices(lattice)) do i
+              return U * interaction_term{i} - mu * N{i}
+              end
     end
 end
 
@@ -343,7 +353,11 @@
 function bose_hubbard_model(elt::Type{<:Number}=ComplexF64,
                             symmetry::Type{<:Sector}=Trivial,
                             lattice::AbstractLattice=InfiniteChain(1);
-                            cutoff::Integer=5, t=1.0, U=1.0, mu=0.0, n::Integer=0)
+                            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)
@@ -352,10 +366,9 @@
     H = @mpoham begin
         sum(nearest_neighbours(lattice)) do (i, j)
             return -t * hopping_term{i,j}
-        end +
-        sum(vertices(lattice)) do i
-            return U / 2 * interaction_term{i} - mu * N{i}
-        end
+        end + sum(vertices(lattice)) do i
+              return U / 2 * interaction_term{i} - mu * N{i}
+              end
     end
 
     if symmetry === Trivial
@@ -370,3 +383,56 @@
 
     return H
 end
+
+#===========================================================================================
+    t-J models
+===========================================================================================#
+
+"""
+    tj_model([elt::Type{<:Number}], [particle_symmetry::Type{<:Sector}],
+                  [spin_symmetry::Type{<:Sector}], [lattice::AbstractLattice];
+                  t, J, mu, sf::Bool=false)
+
+MPO for the hamiltonian of the t-J model, 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`` restrict to the no-double-occupancy subspace. 
+"""
+function tj_model end
+function tj_model(lattice::AbstractLattice; kwargs...)
+    return tj_model(ComplexF64, Trivial, Trivial, lattice; kwargs...)
+end
+function tj_model(particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector};
+                  kwargs...)
+    return tj_model(ComplexF64, particle_symmetry, spin_symmetry; kwargs...)
+end
+function tj_model(elt::Type{<:Number}, lattice::AbstractLattice; kwargs...)
+    return tj_model(elt, Trivial, Trivial, lattice; kwargs...)
+end
+function tj_model(T::Type{<:Number}=ComplexF64,
+                  particle_symmetry::Type{<:Sector}=Trivial,
+                  spin_symmetry::Type{<:Sector}=Trivial,
+                  lattice::AbstractLattice=InfiniteChain(1);
+                  t=2.5,
+                  J=1.0,
+                  mu=0.0,
+                  sf::Bool=false,)
+    hopping = tJ.e_plusmin(T, particle_symmetry, spin_symmetry; sf) +
+              tJ.e_minplus(T, particle_symmetry, spin_symmetry; sf)
+    num = tJ.e_number(T, particle_symmetry, spin_symmetry; sf)
+    heisenberg = tJ.S_exchange(T, particle_symmetry, spin_symmetry; sf) -
+                 (1 / 4) * (num ⊗ num)
+    return @mpoham begin
+        sum(nearest_neighbours(lattice)) do (i, j)
+            return (-t) * hopping{i,j} + J * heisenberg{i,j}
+        end + sum(vertices(lattice)) do i
+            return (-mu) * num{i}
+        end
+    end
+end
+
+# TODO: add (hardcore) bosonic t-J model (https://arxiv.org/abs/2409.15424)

src/models/quantum_chemistry.jl

@@ -42,14 +42,17 @@
                                                          (1, 1 // 2, 1) => 1,
                                                          (2, 0, 0) => 1)
 
-    ap = TensorMap(ones, Elt,
+    ap = TensorMap(ones,
+                   Elt,
                    psp *
                    Vect[(Irrep[U₁] ⊠ Irrep[SU₂] ⊠ FermionParity)]((-1, 1 // 2, 1) => 1),
                    psp)
     blocks(ap)[(U₁(0) ⊠ SU₂(0) ⊠ FermionParity(0))] .*= -sqrt(2)
     blocks(ap)[(U₁(1) ⊠ SU₂(1 // 2) ⊠ FermionParity(1))] .*= 1
 
-    bm = TensorMap(ones, Elt, psp,
+    bm = TensorMap(ones,
+                   Elt,
+                   psp,
                    Vect[(Irrep[U₁] ⊠ Irrep[SU₂] ⊠ FermionParity)]((-1, 1 // 2, 1) => 1) *
                    psp)
     blocks(bm)[(U₁(0) ⊠ SU₂(0) ⊠ FermionParity(0))] .*= sqrt(2)

src/operators/hubbardoperators.jl

@@ -64,7 +64,7 @@ end
 """
     e_plusmin_up(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector})
 
-Return the two-body operator that creates a spin-up electron at the first site and annihilates a spin-up electron at the second.
+Return the two-body operator `e†_{1,↑}, e_{2,↑}` that creates a spin-up electron at the first site and annihilates a spin-up electron at the second.
 """
 e_plusmin_up(P::Type{<:Sector}, S::Type{<:Sector}) = e_plusmin_up(ComplexF64, P, S)
 function e_plusmin_up(T, ::Type{Trivial}, ::Type{Trivial})
@@ -111,7 +111,7 @@ const e⁺e⁻ꜛ = e_plusmin_up
 """
     e_plusmin_down(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector})
 
-Return the two-body operator that creates a spin-down electron at the first site and annihilates a spin-down electron at the second.
+Return the two-body operator `e†_{1,↓}, e_{2,↓}` that creates a spin-down electron at the first site and annihilates a spin-down electron at the second.
 """
 e_plusmin_down(P::Type{<:Sector}, S::Type{<:Sector}) = e_plusmin_down(ComplexF64, P, S)
 function e_plusmin_down(T, ::Type{Trivial}, ::Type{Trivial})
@@ -158,7 +158,9 @@ const e⁺e⁻ꜜ = e_plusmin_down
 """
     e_minplus_up(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector})
 
-Return the two-body operator that annihilates a spin-up electron at the first site and creates a spin-up electron at the second.
+Return the Hermitian conjugate of `e_plusmin_up`, i.e.
+`(e†_{1,↑}, e_{2,↑})† = -e_{1,↑}, e†_{2,↑}` (note the extra minus sign). 
+It annihilates a spin-up electron at the first site and creates a spin-up electron at the second.
 """
 e_minplus_up(P::Type{<:Sector}, S::Type{<:Sector}) = e_minplus_up(ComplexF64, P, S)
 function e_minplus_up(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector})
@@ -169,7 +171,9 @@ const e⁻⁺ꜛ = e_minplus_up
 """
     e_minplus_down(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector})
 
-Return the two-body operator that annihilates a spin-down electron at the first site and creates a spin-down electron at the second.
+Return the Hermitian conjugate of `e_plusmin_down`, i.e.
+`(e†_{1,↓}, e_{2,↓})† = -e_{1,↓}, e†_{2,↓}` (note the extra minus sign). 
+It annihilates a spin-down electron at the first site and creates a spin-down electron at the second.
 """
 e_minplus_down(P::Type{<:Sector}, S::Type{<:Sector}) = e_minplus_down(ComplexF64, P, S)
 function e_minplus_down(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector})