Feature/reduced density matrix (#276)

* 21-july-2024 rDM * 21-july-2024 rDM * 21-july-2024 rDM * 21-july-2024 rDM2 * 21-july-2024 rDM * Update hamiltonians.md * Update hamiltonians.md * Update ReducedDensityMatrix.jl * Update ReducedDensityMatrix.jl * Update ReducedDensityMatrix.jl * Update DensityMatrixDiagonal.jl * Update ReducedDensityMatrix.jl - Changed names of the operators - debugged the code for off-diagonal terms in operators. * Update src/Hamiltonians/ReducedDensityMatrix.jl Co-authored-by: Joachim Brand <joachim.brand@gmail.com> * Update Hamiltonians.jl * Update hamiltonians.md * Update Hamiltonians.jl * Update Hamiltonians.jl * Update DensityMatrixDiagonal.jl * Update Hamiltonians.jl * Update Hamiltonians.jl * Update ReducedDensityMatrix.jl * Update ReducedDensityMatrix.jl * Update ReducedDensityMatrix.jl * Update ReducedDensityMatrix.jl * Update ReducedDensityMatrix.jl * Update src/Hamiltonians/ReducedDensityMatrix.jl Co-authored-by: Joachim Brand <joachim.brand@gmail.com> * Update src/Hamiltonians/ReducedDensityMatrix.jl Co-authored-by: Joachim Brand <joachim.brand@gmail.com> * Update src/Hamiltonians/ReducedDensityMatrix.jl Co-authored-by: Joachim Brand <joachim.brand@gmail.com> * Update src/Hamiltonians/ReducedDensityMatrix.jl Co-authored-by: Joachim Brand <joachim.brand@gmail.com> * Update src/Hamiltonians/ReducedDensityMatrix.jl Co-authored-by: Joachim Brand <joachim.brand@gmail.com> * Update src/Hamiltonians/ReducedDensityMatrix.jl Co-authored-by: Joachim Brand <joachim.brand@gmail.com> * Update src/Hamiltonians/ReducedDensityMatrix.jl Co-authored-by: mtsch <matijacufar@gmail.com> * Update src/Hamiltonians/ReducedDensityMatrix.jl Co-authored-by: Joachim Brand <joachim.brand@gmail.com> * Update ReducedDensityMatrix.jl * Update src/Hamiltonians/Hamiltonians.jl Co-authored-by: Joachim Brand <joachim.brand@gmail.com> * Rename ReducedDensityMatrix.jl to reduced_density_matrix.jl --------- Co-authored-by: Joachim Brand <joachim.brand@gmail.com> Co-authored-by: mtsch <matijacufar@gmail.com>
RimuQMC · Aug 18, 2024 · ecda02c · ecda02c
1 parent c846fb6
commit ecda02c
Show file tree

Hide file tree

Showing 5 changed files with 144 additions and 1 deletion.
diff --git a/docs/src/hamiltonians.md b/docs/src/hamiltonians.md
@@ -80,6 +80,8 @@ G2MomCorrelator
 SuperfluidCorrelator
 StringCorrelator
 DensityMatrixDiagonal
+SingleParticleExcitation
+TwoParticleExcitation
 Momentum
 AxialAngularMomentumHO
 ```

diff --git a/src/Hamiltonians/DensityMatrixDiagonal.jl b/src/Hamiltonians/DensityMatrixDiagonal.jl
@@ -14,6 +14,7 @@ all components is computed.
 
 * [`single_particle_density`](@ref)
 * [`SingleParticleDensity`](@ref)
+* [`SingleParticleExcitation`](@ref)
 """
 struct DensityMatrixDiagonal{C} <: AbstractHamiltonian{Float64}
     mode::Int

diff --git a/src/Hamiltonians/Hamiltonians.jl b/src/Hamiltonians/Hamiltonians.jl
@@ -39,6 +39,8 @@ Other
 - [`G2RealSpace`](@ref)
 - [`G2MomCorrelator`](@ref)
 - [`DensityMatrixDiagonal`](@ref)
+- [`SingleParticleExcitation`](@ref)
+- [`TwoParticleExcitation`](@ref)
 - [`Momentum`](@ref)
 - [`AxialAngularMomentumHO`](@ref)
 
@@ -83,6 +85,7 @@ export FroehlichPolaron
 export ParticleNumberOperator
 
 export G2MomCorrelator, G2RealCorrelator, G2RealSpace, SuperfluidCorrelator, DensityMatrixDiagonal, Momentum
+export SingleParticleExcitation, TwoParticleExcitation
 export StringCorrelator
 
 export CubicGrid, PeriodicBoundaries, HardwallBoundaries, LadderBoundaries
@@ -116,9 +119,9 @@ include("TRSymmetry.jl")
 include("Stoquastic.jl")
 
 include("Transcorrelated1D.jl")
-
 include("correlation_functions.jl")
 include("DensityMatrixDiagonal.jl")
+include("reduced_density_matrix.jl")
 include("Momentum.jl")
 include("particle_number.jl")
 

diff --git a/src/Hamiltonians/reduced_density_matrix.jl b/src/Hamiltonians/reduced_density_matrix.jl
@@ -0,0 +1,104 @@
+"""
+    SingleParticleExcitation(i, j) <: AbstractHamiltonian
+
+Represent the ``{i,j}`` element of the single-particle reduced density matrix:
+
+```math
+\\hat{ρ}^{(1)}_{i,j} = \\hat a^†_{i} \\hat a_{j}
+```
+
+where `i <: Int` and `j <: Int` specify the mode numbers.
+
+# See also
+
+* [`single_particle_density`](@ref)
+* [`SingleParticleDensity`](@ref)
+* [`TwoParticleExcitation`](@ref)
+"""
+struct SingleParticleExcitation{I,J} <: AbstractHamiltonian{Float64}
+end
+
+SingleParticleExcitation(I::Int,J::Int) = SingleParticleExcitation{I,J}()
+
+function Base.show(io::IO, spd::SingleParticleExcitation{I,J}) where {I,J}
+    print(io, "SingleParticleExcitation($(I), $(J))")
+end
+
+LOStructure(::Type{<:SingleParticleExcitation}) = AdjointUnknown()
+
+function diagonal_element(spd::SingleParticleExcitation{I,J}, add::SingleComponentFockAddress) where {I,J}
+    if I != J
+        return 0.0
+    end
+    src = find_mode(add, J)
+    return src.occnum
+end
+
+function num_offdiagonals(spd::SingleParticleExcitation{I,J}, address::SingleComponentFockAddress) where {I,J}
+    if I == J
+        return 0
+    else
+        return 1
+    end
+end
+
+function get_offdiagonal(spd::SingleParticleExcitation{I,J}, add::SingleComponentFockAddress, chosen) where {I,J}
+    src = find_mode(add, J)
+    dst = find_mode(add,I)
+    address, value = excitation(add, (dst,), (src,))
+    return address, value
+end
+
+"""
+    TwoParticleExcitation(i, j, k, l) <: AbstractHamiltonian
+
+Represent the ``{ij, kl}`` element of the two-particle reduced density matrix:
+
+```math
+\\hat{ρ}^{(2)}_{ij, kl} =  \\hat a^†_{i} \\hat a^†_{j} \\hat a_{l} \\hat a_{k} 
+```
+
+where `i`, `j`, `k`, and `l` (all `<: Int`) specify the mode numbers.
+
+# See also
+
+* [`single_particle_density`](@ref)
+* [`SingleParticleDensity`](@ref)
+* [`SingleParticleExcitation`](@ref)
+"""
+struct TwoParticleExcitation{I,J,K,L} <: AbstractHamiltonian{Float64}
+end
+
+TwoParticleExcitation(I::Int,J::Int,K::Int,L::Int) = TwoParticleExcitation{I,J,K,L}()
+
+function Base.show(io::IO, spd::TwoParticleExcitation{I,J,K,L}) where {I,J,K,L}
+    print(io, "TwoParticleExcitation($(I), $(J), $(K), $(L))")
+end
+
+LOStructure(::Type{<:TwoParticleExcitation}) = AdjointUnknown()
+
+function diagonal_element(spd::TwoParticleExcitation{I,J,K,L}, add::SingleComponentFockAddress) where {I,J,K,L}
+    if (I, J) == (K, L) || (I, J) == (L, K)
+        src = find_mode(add, (L, K))
+        dst = find_mode(add,(I, J))
+        address, value = excitation(add, dst, src)
+        return value
+    else
+        return 0.0
+    end
+end
+
+function num_offdiagonals(spd::TwoParticleExcitation{I,J,K,L}, address::SingleComponentFockAddress) where {I,J,K,L}
+    if (I, J) == (K, L)
+        return 0
+    else
+        return 1
+    end
+end
+
+function get_offdiagonal(spd::TwoParticleExcitation{I,J,K,L}, add::SingleComponentFockAddress, chosen) where {I,J,K,L}
+    src = find_mode(add, (L, K))
+    dst = find_mode(add,(I, J))
+    address, value = excitation(add, (dst...,), (src...,))
+    return address, value
+end
diff --git a/test/Hamiltonians.jl b/test/Hamiltonians.jl
@@ -1151,6 +1151,39 @@ using Rimu.Hamiltonians: circshift_dot
         @test LOStructure(DensityMatrixDiagonal(2)) == IsDiagonal()
         @test DensityMatrixDiagonal(15)' === DensityMatrixDiagonal(15)
     end
+
+    @testset "Reduced Density Matrix" begin
+        addr_bose = BoseFS(6,5,4,3,2,1)
+        addr_fermi = FermiFS(1,0,1,0,1,0)
+        for i in 1:6
+            naddr_fermi1, value_f1 = excitation(addr_fermi,find_mode(addr_fermi,(1,)),find_mode(addr_fermi,(i,)))
+            naddr_fermi2, value_f2 = excitation(addr_fermi,find_mode(addr_fermi,(2,)),find_mode(addr_fermi,(i,)))
+            @test diagonal_element(SingleParticleExcitation(i, i), addr_bose) == 7-i
+            @test diagonal_element(SingleParticleExcitation(i, i), addr_fermi) == i%2
+            @test get_offdiagonal(SingleParticleExcitation(2,i),addr_fermi,1) == (naddr_fermi2,value_f2)
+            if i == 1
+                @test diagonal_element(TwoParticleExcitation(1,i,1,i), addr_bose) == 5*(7-i)
+                @test diagonal_element(TwoParticleExcitation(1,i,1,i), addr_fermi) == 0.0
+                @test get_offdiagonal(SingleParticleExcitation(1,i),addr_fermi,1) == (naddr_fermi1,1)
+            else
+                naddr_bose1, value_b1 = excitation(addr_bose,find_mode(addr_bose,(1,)),find_mode(addr_bose,(i,)))
+                naddr_fermi11, value_f11 = excitation(addr_fermi,find_mode(addr_fermi,(1,2,)),find_mode(addr_fermi,(i,1,)))
+                @test diagonal_element(TwoParticleExcitation(1,i,1,i), addr_bose) == 6*(7-i)
+                @test diagonal_element(TwoParticleExcitation(1,i,1,i), addr_fermi) == i%2
+                @test get_offdiagonal(SingleParticleExcitation(1,i),addr_bose,1) == (naddr_bose1,value_b1)
+                @test get_offdiagonal(SingleParticleExcitation(1,i),addr_fermi,1) == (naddr_fermi1,0.0)
+                @test get_offdiagonal(TwoParticleExcitation(1,2,1,i),addr_fermi,1) == (naddr_fermi11,value_f11)
+                if i!=2
+                    naddr_bose11, value_b11 = excitation(addr_bose,find_mode(addr_bose,(1,2,)),find_mode(addr_bose,(i,1,)))
+                    @test get_offdiagonal(TwoParticleExcitation(1,2,1,i),addr_bose,1) == (naddr_bose11,value_b11)
+                end
+            end
+        end
+        @test num_offdiagonals(SingleParticleExcitation(1,2), addr_bose) == 1
+        @test LOStructure(SingleParticleExcitation(1,2)) == AdjointUnknown()
+        @test num_offdiagonals(TwoParticleExcitation(1,2,2,1), addr_bose) == 1
+        @test LOStructure(TwoParticleExcitation(1,2,2,1)) == AdjointUnknown()
+    end
 end
 
 @testset "HubbardReal1DEP" begin