move more code across from Rimu

Project.toml

@@ -7,7 +7,7 @@ version = "0.1.0"
 Rimu = "c53c40cc-bd84-11e9-2cf4-a9fde2b9386e"
 
 [compat]
-Rimu = "0.13.1"
+Rimu = "0.13, 0.14"
 julia = "1.9"
 
 [extras]

src/BoseFS2C.jl

@@ -41,9 +41,48 @@ function Rimu.BitStringAddresses.time_reverse(c::BoseFS2C{NA,NA,M,S,S,N}) where
     return  BoseFS2C{NA,NA,M,S,S,N}(c.bsb, c.bsa)
 end
 
+function Rimu.Hamiltonians.check_tr_address(addr::BoseFS2C{NA,NB,M,SA,SB,N}) where {NA,NB,M,SA,SB,N}
+    if NA ≠ NB || SA ≠ SB
+        throw(ArgumentError("Two component address with equal particle numbers required for `TimeReversalSymmetry`."))
+    end
+end
+
+function Rimu.Interfaces.diagonal_element(dmd::DensityMatrixDiagonal{0}, add::BoseFS2C)
+    return float(find_mode(add.bsa, dmd.mode).occnum + find_mode(add.bsb, dmd.mode).occnum)
+end
+function Rimu.Interfaces.diagonal_element(dmd::DensityMatrixDiagonal{1}, add::BoseFS2C)
+    comp = add.bsa
+    return float(find_mode(comp, dmd.mode).occnum)
+end
+function Rimu.Interfaces.diagonal_element(dmd::DensityMatrixDiagonal{2}, add::BoseFS2C)
+    comp = add.bsb
+    return float(find_mode(comp, dmd.mode).occnum)
+end
+
 function Rimu.Hamiltonians.dimension(::Type{<:BoseFS2C{NA,NB,M}}) where {NA,NB,M}
     return dimension(BoseFS{NA,M}) * dimension(BoseFS{NB,M})
 end
+function Rimu.Hamiltonians.number_conserving_dimension(addr::BoseFS2C)
+    return number_conserving_dimension(addr.bsa) * number_conserving_dimension(addr.bsb)
+end
 
 BoseFS2C(fs::CompositeFS{2}) = BoseFS2C(fs.components...)
 Rimu.CompositeFS(fs::BoseFS2C) = CompositeFS(fs.bsa, fs.bsb)
+
+function Rimu.Interfaces.diagonal_element(m::Momentum{0}, a::BoseFS2C)
+    Rimu.Hamiltonians._momentum((a.bsa, a.bsb), m.fold)
+end
+function Rimu.Interfaces.diagonal_element(m::Momentum{1}, a::BoseFS2C)
+    Rimu.Hamiltonians._momentum(a.bsa, m.fold)
+end
+function Rimu.Interfaces.diagonal_element(m::Momentum{2}, a::BoseFS2C)
+    Rimu.Hamiltonians._momentum(a.bsb, m.fold)
+end
+
+function Rimu.single_particle_density(add::Union{BoseFS2C}; component=0)
+    if component == 0
+        return float.(Tuple(sum(onr(add))))
+    else
+        return float.(Tuple(onr(add)[component]))
+    end
+end

src/G2MomCorrelator.jl

@@ -0,0 +1,160 @@
+
+import Rimu.Hamiltonians: num_offdiagonals, diagonal_element, get_offdiagonal
+
+"""
+    G2MomCorrelator(d::Int,c=:cross) <: AbstractOperator{ComplexF64}
+
+Two-body correlation operator representing the density-density
+correlation at distance `d` of a two component system
+in a momentum-space Fock-state basis.
+It returns a `Complex` value.
+
+Correlation across two components:
+```math
+\\hat{G}^{(2)}(d) = \\frac{1}{M}\\sum_{spqr=1}^M e^{-id(p-q)2π/M} a^†_{s} b^†_{p}  b_q a_r δ_{s+p,q+r}
+```
+Correlation within a single component:
+```math
+\\hat{G}^{(2)}(d) = \\frac{1}{M}\\sum_{spqr=1}^M e^{-id(p-q)2π/M} a^†_{s} a^†_{p}  a_q a_r δ_{s+p,q+r}
+```
+
+The diagonal element, where `(p-q)=0`, is
+```math
+\\frac{1}{M}\\sum_{k,p=1}^M a^†_{k} b^†_{p}  b_p a_k .
+```
+
+# Arguments
+- `d::Integer`: the distance between two particles.
+- `c`: possible instructions: `:cross`: default instruction, computing correlation between
+  particles across two components;
+  `:first`: computing correlation between particles within the first component;
+  `:second`: computing correlation between particles within the second component.
+  These are the only defined instructions, using anything else will produce errors.
+
+# To use on a one-component system
+
+For a system with only one component, e.g. with `BoseFS`, the second argument `c` is
+irrelevant and can be any of the above instructions, one could simply skip this argument
+and let it be the default value.
+
+# See also
+
+* [`BoseHubbardMom1D2C`](@ref)
+* [`BoseFS2C`](@ref)
+* [`Rimu.G2RealCorrelator`](@extref)
+* [`Rimu.G2RealSpace`](@extref)
+* [`Rimu.AbstractOperator`](@extref)
+* [`Rimu.AllOverlaps`](@extref)
+"""
+struct G2MomCorrelator{C} <: AbstractOperator{ComplexF64}
+    d::Int
+
+    function G2MomCorrelator(d, c=:cross)
+        if c == :first
+            return new{1}(d)
+        elseif c == :second
+            return new{2}(d)
+        elseif c == :cross
+            return new{3}(d)
+        else
+            throw(ArgumentError("Unknown instruction for G2MomCorrelator!"))
+        end
+    end
+end
+
+function Rimu.Interfaces.allows_address_type(g2m::G2MomCorrelator, ::Type{A}) where {A}
+    return num_modes(A) > g2m.d
+end
+
+function Base.show(io::IO, g::G2MomCorrelator{C}) where {C}
+    if C == 1
+        print(io, "G2MomCorrelator($(g.d),:first)")
+    elseif C == 2
+        print(io, "G2MomCorrelator($(g.d),:second)")
+    elseif C == 3
+        print(io, "G2MomCorrelator($(g.d),:cross)")
+    end
+end
+
+function num_offdiagonals(g::G2MomCorrelator, addr::SingleComponentFockAddress)
+    m = num_modes(addr)
+    singlies, doublies = num_singly_doubly_occupied_sites(addr)
+    return singlies * (singlies - 1) * (m - 2) + doublies * (m - 1)
+end
+
+num_offdiagonals(g::G2MomCorrelator{1},addr::BoseFS2C) = num_offdiagonals(g, addr.bsa)
+num_offdiagonals(g::G2MomCorrelator{2},addr::BoseFS2C) = num_offdiagonals(g, addr.bsb)
+
+function num_offdiagonals(g::G2MomCorrelator{3}, addr::BoseFS2C)
+    m = num_modes(addr)
+    sa = num_occupied_modes(addr.bsa)
+    sb = num_occupied_modes(addr.bsb)
+    return sa*(m-1)*sb
+end
+
+
+diagonal_element(g::G2MomCorrelator{1}, addr::BoseFS2C) = diagonal_element(g, addr.bsa)
+diagonal_element(g::G2MomCorrelator{2}, addr::BoseFS2C) = diagonal_element(g, addr.bsb)
+
+function diagonal_element(g::G2MomCorrelator{3}, addr::BoseFS2C{NA,NB,M,AA,AB}) where {NA,NB,M,AA,AB}
+    onrep_a = onr(addr.bsa)
+    onrep_b = onr(addr.bsb)
+    gd = 0
+    for p in 1:M
+        iszero(onrep_b[p]) && continue
+        for k in 1:M
+            gd += onrep_a[k] * onrep_b[p] # b†_p b_p a†_k a_k
+        end
+    end
+    return ComplexF64(gd / M)
+end
+
+function diagonal_element(g::G2MomCorrelator, addr::SingleComponentFockAddress)
+    M = num_modes(addr)
+    onrep = onr(addr)
+    gd = 0
+    for p in 1:M
+        iszero(onrep[p]) && continue
+        for k in 1:M
+            gd += onrep[k] * onrep[p] # a†_p a_p a†_k a_k
+        end
+    end
+    return ComplexF64(gd / M)
+end
+
+function get_offdiagonal(g::G2MomCorrelator{1}, addr::A, chosen)::Tuple{A,ComplexF64} where {A<:BoseFS2C}
+    new_bsa, elem = get_offdiagonal(g, addr.bsa, chosen)
+    return A(new_bsa, addr.bsb), elem
+end
+
+function get_offdiagonal(g::G2MomCorrelator{2}, addr::A, chosen)::Tuple{A,ComplexF64} where {A<:BoseFS2C}
+    new_bsb, elem = get_offdiagonal(g, addr.bsb, chosen)
+    return A(addr.bsa, new_bsb), elem
+end
+
+function get_offdiagonal(
+    g::G2MomCorrelator{3},
+    addr::A,
+    chosen,
+    ma=OccupiedModeMap(addr.bsa),
+    mb=OccupiedModeMap(addr.bsb),
+)::Tuple{A,ComplexF64} where {A<:BoseFS2C}
+
+    m = num_modes(addr)
+    new_bsa, new_bsb, gamma, _, _, Δp = momentum_transfer_excitation(
+        addr.bsa, addr.bsb, chosen, ma, mb
+    )
+    gd = exp(-im * g.d * Δp * 2π / m) * gamma
+    return A(new_bsa, new_bsb), ComplexF64(gd / m)
+end
+
+function get_offdiagonal(
+    g::G2MomCorrelator,
+    addr::A,
+    chosen,
+)::Tuple{A,ComplexF64} where {A<:SingleComponentFockAddress}
+    M = num_modes(addr)
+    new_add, gamma, Δp = momentum_transfer_excitation(addr, chosen, OccupiedModeMap(addr))
+    gd = exp(-im * g.d * Δp * 2π / M) * gamma
+    return new_add, ComplexF64(gd / M)
+end

src/RimuLegacyHamiltonians.jl

@@ -1,19 +1,22 @@
 module RimuLegacyHamiltonians
 
-using Rimu: Rimu, BitStringAddresses
-using Rimu.Interfaces: AbstractHamiltonian, IsHermitian
+using Rimu: Rimu, BitStringAddresses, SingleComponentFockAddress
+using Rimu.Interfaces: AbstractHamiltonian, IsHermitian, AbstractOperator
 using Rimu.BitStringAddresses: num_components, onr, near_uniform, time_reverse,
     print_address, AbstractFockAddress, BoseFS, OccupiedModeMap, CompositeFS
 using Rimu.Hamiltonians: Hamiltonians, dimension, HubbardMom1D, starting_address,
     number_conserving_dimension, LOStructure, num_offdiagonals, diagonal_element,
     num_modes, num_occupied_modes, offdiagonals, get_offdiagonal, AbstractOffdiagonals,
-    momentum_transfer_excitation, hopnextneighbour, HubbardReal1D
+    momentum_transfer_excitation, hopnextneighbour, HubbardReal1D, check_tr_address,
+    num_singly_doubly_occupied_sites, DensityMatrixDiagonal, find_mode,
+    number_conserving_dimension, Momentum
 
 export BoseFS2C
-export BoseHubbardMom1D2C, BoseHubbardReal1D2C
+export BoseHubbardMom1D2C, BoseHubbardReal1D2C, G2MomCorrelator
 
 include("BoseFS2C.jl")
 include("BoseHubbardMom1D2C.jl")
 include("BoseHubbardReal1D2C.jl")
+include("G2MomCorrelator.jl")
 
 end