12429834995

Committed 20 Dec 2024 10:31AM UTC coverage: 94.293% (-0.3%) from 94.567%

Build # 12429834995

Build Type

Pull #300

github

Committed by

web-flow

Commit Message

Update src/strategies_and_params/poststepstrategy.jl

Co-authored-by: Joachim Brand <joachim.brand@gmail.com>

Pull Request Pull Request #300: Fix normalisation in `single_particle_density`

Run Details

3 of 3 new or added lines in 1 file covered. (100.0%)

98 existing lines in 23 files now uncovered.

6906 of 7324 relevant lines covered (94.29%)

13553991.69 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

93.47

/src/BitStringAddresses/bosefs.jl

"""
    BoseFS{N,M,S} <: SingleComponentFockAddress

Address type that represents a Fock state of `N` spinless bosons in `M` modes by wrapping a
[`BitString`](@ref), or a [`SortedParticleList`](@ref). Which is wrapped is chosen
automatically based on the properties of the address.

# Constructors

* `BoseFS{[N,M]}(val::Integer...)`: Create `BoseFS{N,M}` from occupation numbers. This is
  type-stable if the number of modes `M` and the number of particles `N` are provided.
  Otherwise, `M` and `N` are inferred from the arguments.

* `BoseFS{[N,M]}(onr)`: Create `BoseFS{N,M}` from occupation number representation, see
  [`onr`](@ref). This is efficient if `N` and `M` are provided, and `onr` is a
  statically-sized collection, such as a `Tuple` or `SVector`.

* `BoseFS{[N,M]}([M, ]pairs...)`: Provide the number of modes `M` and `mode =>
  occupation_number` pairs. If `M` is provided as a type parameter, it should not be
  provided as the first argument.  Useful for creating sparse addresses. `pairs` can be
  multiple arguments or an iterator of pairs.

* `BoseFS{N,M,S}(bs::S)`: Unsafe constructor. Does not check whether the number of
  particles in `bs` is equal to `N`.

* [`@fs_str`](@ref): Addresses are sometimes printed in a compact manner. This
  representation can also be used as a constructor. See the last example below.

# Examples

```jldoctest
julia> BoseFS{6,5}(0, 1, 2, 3, 0)
BoseFS{6,5}(0, 1, 2, 3, 0)

julia> BoseFS([abs(i - 3) ≤ 1 ? i - 1 : 0 for i in 1:5])
BoseFS{6,5}(0, 1, 2, 3, 0)

julia> BoseFS(5, 2 => 1, 3 => 2, 4 => 3)
BoseFS{6,5}(0, 1, 2, 3, 0)

julia> BoseFS{6,5}(i => i - 1 for i in 2:4)
BoseFS{6,5}(0, 1, 2, 3, 0)

julia> fs"|0 1 2 3 0⟩"
BoseFS{6,5}(0, 1, 2, 3, 0)

julia> fs"|b 5: 2 3 3 4 4 4⟩"
BoseFS{6,5}(0, 1, 2, 3, 0)
```

See also: [`SingleComponentFockAddress`](@ref), [`OccupationNumberFS`](@ref),
[`FermiFS`](@ref), [`CompositeFS`](@ref), [`FermiFS2C`](@ref).
"""
struct BoseFS{N,M,S} <: SingleComponentFockAddress{N,M}
    bs::S
end

@inline function BoseFS{N,M,S}(onr::Union{SVector{M},MVector{M},NTuple{M}}) where {N,M,S}
    @boundscheck begin
        sum(onr) == N || throw(ArgumentError(
            "invalid ONR: $N particles expected, $(sum(onr)) given"
        ))
        if S <: BitString
            B = num_bits(S)
            M + N - 1 == B || throw(ArgumentError(
                "invalid ONR: $B-bit BitString does not fit $N particles in $M modes"
            ))
        elseif S <: SortedParticleList
            N == num_particles(S) && M == num_modes(S) || throw(ArgumentError(
                "invalid ONR: $S does not fit $N particles in $M modes"
            ))
        end
    end
    return BoseFS{N,M,S}(from_bose_onr(S, onr))
end
function BoseFS{N,M}(onr::Union{AbstractArray{<:Integer},NTuple{M,<:Integer}}) where {N,M}
    @boundscheck begin
        sum(onr) == N || throw(ArgumentError(
            "invalid ONR: $N particles expected, $(sum(onr)) given"
        ))
        length(onr) == M || throw(ArgumentError(
            "invalid ONR: $M modes expected, $(length(onr)) given"
        ))
    end
    spl_type = select_int_type(M)

    # Pick smaller address type, but prefer sparse.
    # Alway pick dense if it fits into one chunk.

    # Compute the size of container in words
    sparse_sizeof = ceil(Int, N * sizeof(spl_type) / 8)
    dense_sizeof = ceil(Int, (N + M - 1) / 64)
    if dense_sizeof == 1 || dense_sizeof < sparse_sizeof
        S = typeof(BitString{M + N - 1}(0))
    else
        S = SortedParticleList{N,M,spl_type}
    end
    return BoseFS{N,M,S}(from_bose_onr(S, onr))
end
function BoseFS(onr::Union{AbstractArray,Tuple})
    M = length(onr)
    N = sum(onr)
    return BoseFS{N,M}(onr)
end
BoseFS(vals::Integer...) = BoseFS(vals) # specify occupation numbers
BoseFS(val::Integer) = BoseFS((val,)) # single mode address
BoseFS{N,M}(vals::Integer...) where {N,M} = BoseFS{N,M}(vals)

BoseFS(M::Integer, pairs::Pair...) = BoseFS(M, pairs)
BoseFS(M::Integer, pairs) = BoseFS(sparse_to_onr(M, pairs))
BoseFS{N,M}(pairs::Pair...) where {N,M} = BoseFS{N,M}(pairs)
BoseFS{N,M}(pairs) where {N,M} = BoseFS{N,M}(sparse_to_onr(M, pairs))

function print_address(io::IO, b::BoseFS{N,M}; compact=false) where {N,M}
    if compact && b.bs isa SortedParticleList
        print(io, "|b ", M, ": ", join(Int.(b.bs.storage), ' '), "⟩")
    elseif compact
        print(io, "|", join(onr(b), ' '), "⟩")
    elseif b.bs isa SortedParticleList
        print(io, "BoseFS{$N,$M}(", onr_sparse_string(onr(b)), ")")
    else
        print(io, "BoseFS{$N,$M}", tuple(onr(b)...))
    end
end

Base.bitstring(b::BoseFS) = bitstring(b.bs) # TODO rename?

Base.isless(a::BoseFS, b::BoseFS) = isless(a.bs, b.bs)
Base.hash(bba::BoseFS,  h::UInt) = hash(bba.bs, h)
Base.:(==)(a::BoseFS, b::BoseFS) = a.bs == b.bs

"""
    near_uniform_onr(N, M) -> onr::SVector{M,Int}

Create occupation number representation `onr` distributing `N` particles in `M`
modes in a close-to-uniform fashion with each mode filled with at least
`N ÷ M` particles and at most with `N ÷ M + 1` particles.
"""
function near_uniform_onr(n::Number, m::Number)
    return near_uniform_onr(Val(n),Val(m))
end
function near_uniform_onr(::Val{N}, ::Val{M}) where {N, M}
    fillingfactor, extras = divrem(N, M)
    # startonr = fill(fillingfactor,M)
    startonr = fillingfactor * @MVector ones(Int,M)
    startonr[1:extras] .+= 1
    return SVector{M}(startonr)
end

"""
    near_uniform(BoseFS{N,M}) -> BoseFS{N,M}

Create bosonic Fock state with near uniform occupation number of `M` modes with
a total of `N` particles.

# Examples
```jldoctest
julia> near_uniform(BoseFS{7,5})
BoseFS{7,5}(2, 2, 1, 1, 1)

julia> near_uniform(FermiFS{3,5})
FermiFS{3,5}(1, 1, 1, 0, 0)
```
"""
function near_uniform(::Type{<:BoseFS{N,M}}) where {N,M}
    return BoseFS{N,M}(near_uniform_onr(Val(N),Val(M)))
end
near_uniform(b::AbstractFockAddress) = near_uniform(typeof(b))

onr(b::BoseFS{<:Any,M}) where {M} = to_bose_onr(b.bs, Val(M))
const occupation_number_representation = onr # resides here because `onr` has to be defined

function Base.reverse(b::BoseFS)
    return typeof(b)(reverse(b.bs))
end

# For vacuum state
function num_occupied_modes(b::BoseFS{0})
    return 0
end
function num_occupied_modes(b::BoseFS)
    return bose_num_occupied_modes(b.bs)
end
function occupied_modes(b::BoseFS{N,M,S}) where {N,M,S}
    return BoseOccupiedModes{N,M,S}(b.bs)
end

function find_mode(b::BoseFS, index)
    last_occnum = last_mode = last_offset = 0
    for (occnum, mode, offset) in occupied_modes(b)
        dist = index - mode
        if dist == 0
            return BoseFSIndex(occnum, index, offset)
        elseif dist < 0
            return BoseFSIndex(0, index, offset + dist)
        end
        last_occnum = occnum
        last_mode = mode
        last_offset = offset
    end
    offset = last_offset + last_occnum + index - last_mode
    return BoseFSIndex(0, index, offset)
end
# Multiple in a single pass
function find_mode(b::BoseFS, indices::NTuple{N}) where {N}
    # Idea: find permutation, then use the permutation to find indices in order even though
    # they are not sorted.
    perm = sortperm(SVector(indices))
    # perm_i is the index in permutation and goes from 1:N.
    perm_i = 1
    # curr_i points to indices and result
    curr_i = perm[1]
    # index is the current index we are looking for.
    index = indices[curr_i]

    result = ntuple(_ -> BoseFSIndex(0, 0, 0), Val(N))
    last_occnum = last_mode = last_offset = 0
    @inbounds for (occnum, mode, offset) in occupied_modes(b)
        dist = index - mode
        # While loop handles duplicate entries in indices.
        while dist ≤ 0
            if dist == 0
                @set! result[curr_i] = BoseFSIndex(occnum, mode, offset)
            else
                @set! result[curr_i] = BoseFSIndex(0, index, offset + dist)
            end
            perm_i += 1
            perm_i > N && return result
            curr_i = perm[perm_i]
            index = indices[curr_i]
            dist = index - mode
        end
        last_occnum = occnum
        last_mode = mode
        last_offset = offset
    end
    # Now we have to find all indices that appear after the last occupied site.
    # While true because we break out of the loop early anyway.
    @inbounds while true
        offset = last_offset + last_occnum + index - last_mode
        @set! result[curr_i] = BoseFSIndex(0, index, offset)
        perm_i += 1
        perm_i > N && return result
        curr_i = perm[perm_i]
        index = indices[curr_i]
    end
    return result # not reached
end

# find_occupied_mode provided by generic implementation

function excitation(b::B, creations, destructions) where {B<:BoseFS}
    new_bs, val = bose_excitation(b.bs, creations, destructions)
    return B(new_bs), val
end

"""
    new_address, value = hopnextneighbour(add, chosen, boundary_condition)

Compute the new address of a hopping event for the Hubbard model. Returns the new
address and the square root of product of occupation numbers of the involved modes
multiplied by a term consistent with boundary condition as the `value`. 
The following boundary conditions are supported:

* `:periodic`: hopping over the boundary gives does not change the `value`.
* `:twisted`: hopping over the boundary flips the sign of the `value`.
* `:hard_wall`: hopping over the boundary gives a `value` of zero.
* `θ <: Number`: hopping over the boundary gives a `value` multiplied by ``\\exp(iθ)`` or ``\\exp(−iθ)`` depending on the direction of hopping.

The off-diagonals are indexed as follows:

* `(chosen + 1) ÷ 2` selects the hopping site.
* Even `chosen` indicates a hop to the left.
* Odd `chosen` indicates a hop to the right.

# Example

```jldoctest
julia> using Rimu.Hamiltonians: hopnextneighbour

julia> hopnextneighbour(BoseFS(1, 0, 1), 3)
(BoseFS{2,3}(2, 0, 0), 1.4142135623730951)

julia> hopnextneighbour(BoseFS(1, 0, 1), 4)
(BoseFS{2,3}(1, 1, 0), 1.0)

julia> hopnextneighbour(BoseFS(1, 0, 1), 3, :twisted)
(BoseFS{2,3}(2, 0, 0), -1.4142135623730951)

julia> hopnextneighbour(BoseFS(1, 0, 1), 3, :hard_wall)
(BoseFS{2,3}(2, 0, 0), 0.0)

julia> hopnextneighbour(BoseFS(1, 0, 1), 3, π/4)
(BoseFS{2,3}(2, 0, 0), 1.0000000000000002 + 1.0im)
```
"""
function hopnextneighbour(b::BoseFS{N,M,A}, chosen) where {N,M,A<:BitString}
    address = b.bs
    T = chunk_type(address)
    site = (chosen + 1) >>> 0x1
    if isodd(chosen) # Hopping to the right
        next = 0
        curr = 0
        offset = 0
        sc = 0
        reached_end = false
        for (i, (num, sn, bit)) in enumerate(occupied_modes(b))
            next = num * (sn == sc + 1) # only set next to > 0 if sites are neighbours
            reached_end = i == site + 1
            reached_end && break
            curr = num
            offset = bit + num
            sc = sn
        end
        if sc == M
            new_address = (address << 0x1) | A(T(1))
            prod = curr * (trailing_ones(address) + 1) # mul occupation num of first obital
        else
            next *= reached_end
            new_address = address ⊻ A(T(3)) << ((offset - 1) % T)
            prod = curr * (next + 1)
        end
    else # Hopping to the left
        if site == 1 && isodd(address)
            # For leftmost site, we shift the whole address circularly by one bit.
            new_address = (address >>> 0x1) | A(T(1)) << ((N + M - 2) % T)
            prod = trailing_ones(address) * leading_ones(new_address)
        else
            prev = 0
            curr = 0
            offset = 0
            sp = 0
            for (i, (num, sc, bit)) in enumerate(occupied_modes(b))
                prev = curr * (sc == sp + 1) # only set prev to > 0 if sites are neighbours
                curr = num
                offset = bit
                i == site && break
                sp = sc
            end
            new_address = address ⊻ A(T(3)) << ((offset - 1) % T)
            prod = curr * (prev + 1)
        end
    end
    return BoseFS{N,M,A}(new_address), √prod
end

function hopnextneighbour(b::SingleComponentFockAddress, i)
    src = find_occupied_mode(b, (i + 1) >>> 0x1)
    dst = find_mode(b, mod1(src.mode + ifelse(isodd(i), 1, -1), num_modes(b)))

    new_b, val = excitation(b, (dst,), (src,))
    return new_b, val
end

function hopnextneighbour(
    b::SingleComponentFockAddress, i, boundary_condition::Symbol)
    src = find_occupied_mode(b, (i + 1) >>> 0x1)
    dir = ifelse(isodd(i), 1, -1)
    dst = find_mode(b, mod1(src.mode + dir, num_modes(b)))
    new_b, val = excitation(b, (dst,), (src,))
    on_boundary = src.mode == 1 && dir == -1 || src.mode == num_modes(b) && dir == 1
    if boundary_condition == :twisted && on_boundary
        return new_b, -val
    elseif boundary_condition == :hard_wall && on_boundary
        return new_b, 0.0
    else
        return new_b, val
    end
end

function hopnextneighbour(b::SingleComponentFockAddress, i, boundary_condition::Number)
    src = find_occupied_mode(b, (i + 1) >>> 0x1)
    dir = ifelse(isodd(i), 1, -1)
    dst = find_mode(b, mod1(src.mode + dir, num_modes(b)))
    new_b, val = excitation(b, (dst,), (src,))
    if (src.mode == 1 && dir == -1)
        return new_b, val*exp(-im*boundary_condition)
    elseif (src.mode == num_modes(b) && dir == 1)
        return new_b, val*exp(im*boundary_condition)
    else
        return new_b, complex(val)
    end
end

"""
    bose_hubbard_interaction(address)

Return ``Σ_i n_i (n_i-1)`` for computing the Bose-Hubbard on-site interaction (without the
``U`` prefactor.)

# Example

```jldoctest
julia> Hamiltonians.bose_hubbard_interaction(BoseFS{4,4}((2,1,1,0)))
2
julia> Hamiltonians.bose_hubbard_interaction(BoseFS{4,4}((3,0,1,0)))
6
```
"""
function bose_hubbard_interaction(b::BoseFS{<:Any,<:Any,A}) where {A<:BitString}
    return bose_hubbard_interaction(Val(num_chunks(A)), b)
end
function bose_hubbard_interaction(b::SingleComponentFockAddress)
    return bose_hubbard_interaction(nothing, b)
end

@inline function bose_hubbard_interaction(_, b::SingleComponentFockAddress)
    result = 0
    for (n, _, _) in occupied_modes(b)
        result += n * (n - 1)
    end
    return result
end

@inline function bose_hubbard_interaction(::Val{1}, b::BoseFS{<:Any,<:Any,<:BitString})
    # currently this ammounts to counting occupation numbers of modes
    chunk = chunks(b.bs)[1]
    matrixelementint = 0
    while !iszero(chunk)
        chunk >>>= (trailing_zeros(chunk) % UInt) # proceed to next occupied mode
        bosonnumber = trailing_ones(chunk) # count how many bosons inside
        # surpsingly it is faster to not check whether this is nonzero and do the
        # following operations anyway
        chunk >>>= (bosonnumber % UInt) # remove the counted mode
        matrixelementint += bosonnumber * (bosonnumber - 1)
    end
    return matrixelementint
end

1	"""
2	BoseFS{N,M,S} <: SingleComponentFockAddress
3
4	Address type that represents a Fock state of `N` spinless bosons in `M` modes by wrapping a
5	[`BitString`](@ref), or a [`SortedParticleList`](@ref). Which is wrapped is chosen
6	automatically based on the properties of the address.
7
8	# Constructors
9
10	* `BoseFS{[N,M]}(val::Integer...)`: Create `BoseFS{N,M}` from occupation numbers. This is
11	type-stable if the number of modes `M` and the number of particles `N` are provided.
12	Otherwise, `M` and `N` are inferred from the arguments.
13
14	* `BoseFS{[N,M]}(onr)`: Create `BoseFS{N,M}` from occupation number representation, see
15	[`onr`](@ref). This is efficient if `N` and `M` are provided, and `onr` is a
16	statically-sized collection, such as a `Tuple` or `SVector`.
17
18	* `BoseFS{[N,M]}([M, ]pairs...)`: Provide the number of modes `M` and `mode =>
19	occupation_number` pairs. If `M` is provided as a type parameter, it should not be
20	provided as the first argument. Useful for creating sparse addresses. `pairs` can be
21	multiple arguments or an iterator of pairs.
22
23	* `BoseFS{N,M,S}(bs::S)`: Unsafe constructor. Does not check whether the number of
24	particles in `bs` is equal to `N`.
25
26	* [`@fs_str`](@ref): Addresses are sometimes printed in a compact manner. This
27	representation can also be used as a constructor. See the last example below.
28
29	# Examples
30
31	```jldoctest
32	julia> BoseFS{6,5}(0, 1, 2, 3, 0)
33	BoseFS{6,5}(0, 1, 2, 3, 0)
34
35	julia> BoseFS([abs(i - 3) ≤ 1 ? i - 1 : 0 for i in 1:5])
36	BoseFS{6,5}(0, 1, 2, 3, 0)
37
38	julia> BoseFS(5, 2 => 1, 3 => 2, 4 => 3)
39	BoseFS{6,5}(0, 1, 2, 3, 0)
40
41	julia> BoseFS{6,5}(i => i - 1 for i in 2:4)
42	BoseFS{6,5}(0, 1, 2, 3, 0)
43
44	julia> fs"\|0 1 2 3 0⟩"
45	BoseFS{6,5}(0, 1, 2, 3, 0)
46
47	julia> fs"\|b 5: 2 3 3 4 4 4⟩"
48	BoseFS{6,5}(0, 1, 2, 3, 0)
49	```
50
51	See also: [`SingleComponentFockAddress`](@ref), [`OccupationNumberFS`](@ref),
52	[`FermiFS`](@ref), [`CompositeFS`](@ref), [`FermiFS2C`](@ref).
53	"""
54	struct BoseFS{N,M,S} <: SingleComponentFockAddress{N,M}
55	bs::S	287,007,474✔
56	end
57
UNCOV 58	@inline function BoseFS{N,M,S}(onr::Union{SVector{M},MVector{M},NTuple{M}}) where {N,M,S}	×
UNCOV 59	@boundscheck begin	×
UNCOV 60	sum(onr) == N \|\| throw(ArgumentError(	×
61	"invalid ONR: $N particles expected, $(sum(onr)) given"
62	))
UNCOV 63	if S <: BitString	×
UNCOV 64	B = num_bits(S)	×
UNCOV 65	M + N - 1 == B \|\| throw(ArgumentError(	×
66	"invalid ONR: $B-bit BitString does not fit $N particles in $M modes"
67	))
68	elseif S <: SortedParticleList	×
69	N == num_particles(S) && M == num_modes(S) \|\| throw(ArgumentError(	×
70	"invalid ONR: $S does not fit $N particles in $M modes"
71	))
72	end
73	end
UNCOV 74	return BoseFS{N,M,S}(from_bose_onr(S, onr))	×
75	end
76	function BoseFS{N,M}(onr::Union{AbstractArray{<:Integer},NTuple{M,<:Integer}}) where {N,M}	544,228✔
77	@boundscheck begin	544,228✔
78	sum(onr) == N \|\| throw(ArgumentError(	544,229✔
79	"invalid ONR: $N particles expected, $(sum(onr)) given"
80	))
81	length(onr) == M \|\| throw(ArgumentError(	544,228✔
82	"invalid ONR: $M modes expected, $(length(onr)) given"
83	))
84	end
85	spl_type = select_int_type(M)	544,226✔
86
87	# Pick smaller address type, but prefer sparse.
88	# Alway pick dense if it fits into one chunk.
89
90	# Compute the size of container in words
91	sparse_sizeof = ceil(Int, N * sizeof(spl_type) / 8)	544,226✔
92	dense_sizeof = ceil(Int, (N + M - 1) / 64)	544,226✔
93	if dense_sizeof == 1 \|\| dense_sizeof < sparse_sizeof	544,226✔
94	S = typeof(BitString{M + N - 1}(0))	539,458✔
95	else
96	S = SortedParticleList{N,M,spl_type}	4,768✔
97	end
98	return BoseFS{N,M,S}(from_bose_onr(S, onr))	544,951✔
99	end
100	function BoseFS(onr::Union{AbstractArray,Tuple})	309,022✔
101	M = length(onr)	309,022✔
102	N = sum(onr)	309,882✔
103	return BoseFS{N,M}(onr)	309,731✔
104	end
105	BoseFS(vals::Integer...) = BoseFS(vals) # specify occupation numbers	116✔
106	BoseFS(val::Integer) = BoseFS((val,)) # single mode address	2✔
107	BoseFS{N,M}(vals::Integer...) where {N,M} = BoseFS{N,M}(vals)	80✔
108
109	BoseFS(M::Integer, pairs::Pair...) = BoseFS(M, pairs)	34✔
110	BoseFS(M::Integer, pairs) = BoseFS(sparse_to_onr(M, pairs))	562✔
111	BoseFS{N,M}(pairs::Pair...) where {N,M} = BoseFS{N,M}(pairs)	11✔
112	BoseFS{N,M}(pairs) where {N,M} = BoseFS{N,M}(sparse_to_onr(M, pairs))	12✔
113
114	function print_address(io::IO, b::BoseFS{N,M}; compact=false) where {N,M}	1,802✔
115	if compact && b.bs isa SortedParticleList	901✔
116	print(io, "\|b ", M, ": ", join(Int.(b.bs.storage), ' '), "⟩")	10✔
117	elseif compact	891✔
118	print(io, "\|", join(onr(b), ' '), "⟩")	546✔
119	elseif b.bs isa SortedParticleList	345✔
120	print(io, "BoseFS{$N,$M}(", onr_sparse_string(onr(b)), ")")	40✔
121	else
122	print(io, "BoseFS{$N,$M}", tuple(onr(b)...))	305✔
123	end
124	end
125
126	Base.bitstring(b::BoseFS) = bitstring(b.bs) # TODO rename?	×
127
128	Base.isless(a::BoseFS, b::BoseFS) = isless(a.bs, b.bs)	136,543✔
129	Base.hash(bba::BoseFS, h::UInt) = hash(bba.bs, h)	763,107,614✔
130	Base.:(==)(a::BoseFS, b::BoseFS) = a.bs == b.bs	49,503,401✔
131
132	"""
133	near_uniform_onr(N, M) -> onr::SVector{M,Int}
134
135	Create occupation number representation `onr` distributing `N` particles in `M`
136	modes in a close-to-uniform fashion with each mode filled with at least
137	`N ÷ M` particles and at most with `N ÷ M + 1` particles.
138	"""
139	function near_uniform_onr(n::Number, m::Number)	×
140	return near_uniform_onr(Val(n),Val(m))	×
141	end
142	function near_uniform_onr(::Val{N}, ::Val{M}) where {N, M}	21✔
143	fillingfactor, extras = divrem(N, M)	21✔
144	# startonr = fill(fillingfactor,M)
145	startonr = fillingfactor * @MVector ones(Int,M)	21✔
146	startonr[1:extras] .+= 1	21✔
147	return SVector{M}(startonr)	21✔
148	end
149
150	"""
151	near_uniform(BoseFS{N,M}) -> BoseFS{N,M}
152
153	Create bosonic Fock state with near uniform occupation number of `M` modes with
154	a total of `N` particles.
155
156	# Examples
157	```jldoctest
158	julia> near_uniform(BoseFS{7,5})
159	BoseFS{7,5}(2, 2, 1, 1, 1)
160
161	julia> near_uniform(FermiFS{3,5})
162	FermiFS{3,5}(1, 1, 1, 0, 0)
163	```
164	"""
165	function near_uniform(::Type{<:BoseFS{N,M}}) where {N,M}	21✔
166	return BoseFS{N,M}(near_uniform_onr(Val(N),Val(M)))	21✔
167	end
168	near_uniform(b::AbstractFockAddress) = near_uniform(typeof(b))	1✔
169
170	onr(b::BoseFS{<:Any,M}) where {M} = to_bose_onr(b.bs, Val(M))	32,710,437✔
171	const occupation_number_representation = onr # resides here because `onr` has to be defined
172
173	function Base.reverse(b::BoseFS)	1,766✔
174	return typeof(b)(reverse(b.bs))	1,766✔
175	end
176
177	# For vacuum state
178	function num_occupied_modes(b::BoseFS{0})	1✔
179	return 0	1✔
180	end
181	function num_occupied_modes(b::BoseFS)	16,716,046✔
182	return bose_num_occupied_modes(b.bs)	68,225,440✔
183	end
184	function occupied_modes(b::BoseFS{N,M,S}) where {N,M,S}	499,450,754✔
185	return BoseOccupiedModes{N,M,S}(b.bs)	499,455,629✔
186	end
187
188	function find_mode(b::BoseFS, index)	11,115,559✔
189	last_occnum = last_mode = last_offset = 0	11,115,559✔
190	for (occnum, mode, offset) in occupied_modes(b)	20,399,812✔
191	dist = index - mode	64,979,921✔
192	if dist == 0	64,979,921✔
193	return BoseFSIndex(occnum, index, offset)	5,900,599✔
194	elseif dist < 0	59,079,322✔
195	return BoseFSIndex(0, index, offset + dist)	3,777,942✔
196	end
197	last_occnum = occnum	55,301,380✔
198	last_mode = mode	55,301,380✔
199	last_offset = offset	55,301,380✔
200	end	109,170,789✔
201	offset = last_offset + last_occnum + index - last_mode	1,437,018✔
202	return BoseFSIndex(0, index, offset)	1,437,018✔
203	end
204	# Multiple in a single pass
205	function find_mode(b::BoseFS, indices::NTuple{N}) where {N}	228,787,346✔
206	# Idea: find permutation, then use the permutation to find indices in order even though
207	# they are not sorted.
208	perm = sortperm(SVector(indices))	228,787,358✔
209	# perm_i is the index in permutation and goes from 1:N.
210	perm_i = 1	228,787,363✔
211	# curr_i points to indices and result
212	curr_i = perm[1]	228,787,355✔
213	# index is the current index we are looking for.
214	index = indices[curr_i]	228,787,362✔
215
216	result = ntuple(_ -> BoseFSIndex(0, 0, 0), Val(N))	686,361,908✔
217	last_occnum = last_mode = last_offset = 0	228,787,426✔
218	@inbounds for (occnum, mode, offset) in occupied_modes(b)	457,574,777✔
219	dist = index - mode	640,722,782✔
220	# While loop handles duplicate entries in indices.
221	while dist ≤ 0	855,250,525✔
222	if dist == 0	380,016,127✔
223	@set! result[curr_i] = BoseFSIndex(occnum, mode, offset)	193,298,865✔
224	else
225	@set! result[curr_i] = BoseFSIndex(0, index, offset + dist)	186,717,380✔
226	end
227	perm_i += 1	380,016,187✔
228	perm_i > N && return result	380,016,197✔
229	curr_i = perm[perm_i]	214,528,314✔
230	index = indices[curr_i]	214,528,313✔
231	dist = index - mode	214,528,319✔
232	end	214,528,330✔
233	last_occnum = occnum	475,235,157✔
234	last_mode = mode	475,235,185✔
235	last_offset = offset	475,235,208✔
236	end	887,171,038✔
237	# Now we have to find all indices that appear after the last occupied site.
238	# While true because we break out of the loop early anyway.
239	@inbounds while true	63,299,374✔
240	offset = last_offset + last_occnum + index - last_mode	77,558,436✔
241	@set! result[curr_i] = BoseFSIndex(0, index, offset)	77,558,438✔
242	perm_i += 1	77,558,441✔
243	perm_i > N && return result	77,558,442✔
244	curr_i = perm[perm_i]	14,259,063✔
245	index = indices[curr_i]	14,259,063✔
246	end	14,259,063✔
247	return result # not reached	×
248	end
249
250	# find_occupied_mode provided by generic implementation
251
252	function excitation(b::B, creations, destructions) where {B<:BoseFS}	231,721,474✔
253	new_bs, val = bose_excitation(b.bs, creations, destructions)	234,095,956✔
254	return B(new_bs), val	231,721,526✔
255	end
256
257	"""
258	new_address, value = hopnextneighbour(add, chosen, boundary_condition)
259
260	Compute the new address of a hopping event for the Hubbard model. Returns the new
261	address and the square root of product of occupation numbers of the involved modes
262	multiplied by a term consistent with boundary condition as the `value`.
263	The following boundary conditions are supported:
264
265	* `:periodic`: hopping over the boundary gives does not change the `value`.
266	* `:twisted`: hopping over the boundary flips the sign of the `value`.
267	* `:hard_wall`: hopping over the boundary gives a `value` of zero.
268	* `θ <: Number`: hopping over the boundary gives a `value` multiplied by ``\\exp(iθ)`` or ``\\exp(−iθ)`` depending on the direction of hopping.
269
270	The off-diagonals are indexed as follows:
271
272	* `(chosen + 1) ÷ 2` selects the hopping site.
273	* Even `chosen` indicates a hop to the left.
274	* Odd `chosen` indicates a hop to the right.
275
276	# Example
277
278	```jldoctest
279	julia> using Rimu.Hamiltonians: hopnextneighbour
280
281	julia> hopnextneighbour(BoseFS(1, 0, 1), 3)
282	(BoseFS{2,3}(2, 0, 0), 1.4142135623730951)
283
284	julia> hopnextneighbour(BoseFS(1, 0, 1), 4)
285	(BoseFS{2,3}(1, 1, 0), 1.0)
286
287	julia> hopnextneighbour(BoseFS(1, 0, 1), 3, :twisted)
288	(BoseFS{2,3}(2, 0, 0), -1.4142135623730951)
289
290	julia> hopnextneighbour(BoseFS(1, 0, 1), 3, :hard_wall)
291	(BoseFS{2,3}(2, 0, 0), 0.0)
292
293	julia> hopnextneighbour(BoseFS(1, 0, 1), 3, π/4)
294	(BoseFS{2,3}(2, 0, 0), 1.0000000000000002 + 1.0im)
295	```
296	"""
297	function hopnextneighbour(b::BoseFS{N,M,A}, chosen) where {N,M,A<:BitString}	54,619,717✔
298	address = b.bs	54,645,938✔
299	T = chunk_type(address)	54,664,971✔
300	site = (chosen + 1) >>> 0x1	54,686,076✔
301	if isodd(chosen) # Hopping to the right	54,688,909✔
302	next = 0	27,400,211✔
303	curr = 0	27,402,068✔
304	offset = 0	27,403,276✔
305	sc = 0	27,404,157✔
306	reached_end = false	27,402,220✔
307	for (i, (num, sn, bit)) in enumerate(occupied_modes(b))	54,600,746✔
308	next = num * (sn == sc + 1) # only set next to > 0 if sites are neighbours	91,079,366✔
309	reached_end = i == site + 1	91,054,164✔
310	reached_end && break	91,088,139✔
311	curr = num	71,028,251✔
312	offset = bit + num	71,037,834✔
313	sc = sn	71,049,279✔
314	end	71,043,616✔
315	if sc == M	27,348,852✔
316	new_address = (address << 0x1) \| A(T(1))	3,370,805✔
317	prod = curr * (trailing_ones(address) + 1) # mul occupation num of first obital	3,371,768✔
318	else
319	next *= reached_end	24,005,551✔
320	new_address = address ⊻ A(T(3)) << ((offset - 1) % T)	24,010,312✔
321	prod = curr * (next + 1)	24,061,418✔
322	end
323	else # Hopping to the left
324	if site == 1 && isodd(address)	27,401,705✔
325	# For leftmost site, we shift the whole address circularly by one bit.
326	new_address = (address >>> 0x1) \| A(T(1)) << ((N + M - 2) % T)	4,266,394✔
327	prod = trailing_ones(address) * leading_ones(new_address)	4,267,743✔
328	else
329	prev = 0	23,149,318✔
330	curr = 0	23,141,679✔
331	offset = 0	23,152,085✔
332	sp = 0	23,152,682✔
333	for (i, (num, sc, bit)) in enumerate(occupied_modes(b))	46,217,203✔
334	prev = curr * (sc == sp + 1) # only set prev to > 0 if sites are neighbours	66,909,553✔
335	curr = num	66,905,384✔
336	offset = bit	66,908,832✔
337	i == site && break	66,896,912✔
338	sp = sc	43,891,099✔
339	end	43,892,701✔
340	new_address = address ⊻ A(T(3)) << ((offset - 1) % T)	23,136,140✔
341	prod = curr * (prev + 1)	23,174,848✔
342	end
343	end
344	return BoseFS{N,M,A}(new_address), √prod	54,734,545✔
345	end
346
347	function hopnextneighbour(b::SingleComponentFockAddress, i)	854✔
348	src = find_occupied_mode(b, (i + 1) >>> 0x1)	854✔
349	dst = find_mode(b, mod1(src.mode + ifelse(isodd(i), 1, -1), num_modes(b)))	1,054✔
350
351	new_b, val = excitation(b, (dst,), (src,))	1,442✔
352	return new_b, val	854✔
353	end
354
355	function hopnextneighbour(	861✔
356	b::SingleComponentFockAddress, i, boundary_condition::Symbol)
357	src = find_occupied_mode(b, (i + 1) >>> 0x1)	2,587✔
358	dir = ifelse(isodd(i), 1, -1)	861✔
359	dst = find_mode(b, mod1(src.mode + dir, num_modes(b)))	2,740✔
360	new_b, val = excitation(b, (dst,), (src,))	1,528✔
361	on_boundary = src.mode == 1 && dir == -1 \|\| src.mode == num_modes(b) && dir == 1	1,596✔
362	if boundary_condition == :twisted && on_boundary	861✔
363	return new_b, -val	13✔
364	elseif boundary_condition == :hard_wall && on_boundary	848✔
365	return new_b, 0.0	9✔
366	else
367	return new_b, val	839✔
368	end
369	end
370
371	function hopnextneighbour(b::SingleComponentFockAddress, i, boundary_condition::Number)	29✔
372	src = find_occupied_mode(b, (i + 1) >>> 0x1)	42✔
373	dir = ifelse(isodd(i), 1, -1)	29✔
374	dst = find_mode(b, mod1(src.mode + dir, num_modes(b)))	29✔
375	new_b, val = excitation(b, (dst,), (src,))	58✔
376	if (src.mode == 1 && dir == -1)	29✔
377	return new_b, valexp(-imboundary_condition)	6✔
378	elseif (src.mode == num_modes(b) && dir == 1)	23✔
379	return new_b, valexp(imboundary_condition)	5✔
380	else
381	return new_b, complex(val)	18✔
382	end
383	end
384
385	"""
386	bose_hubbard_interaction(address)
387
388	Return ``Σ_i n_i (n_i-1)`` for computing the Bose-Hubbard on-site interaction (without the
389	``U`` prefactor.)
390
391	# Example
392
393	```jldoctest
394	julia> Hamiltonians.bose_hubbard_interaction(BoseFS{4,4}((2,1,1,0)))
395	2
396	julia> Hamiltonians.bose_hubbard_interaction(BoseFS{4,4}((3,0,1,0)))
397	6
398	```
399	"""
400	function bose_hubbard_interaction(b::BoseFS{<:Any,<:Any,A}) where {A<:BitString}	47,098,533✔
401	return bose_hubbard_interaction(Val(num_chunks(A)), b)	47,109,571✔
402	end
403	function bose_hubbard_interaction(b::SingleComponentFockAddress)	70✔
404	return bose_hubbard_interaction(nothing, b)	70✔
405	end
406
407	@inline function bose_hubbard_interaction(_, b::SingleComponentFockAddress)	103✔
408	result = 0	103✔
409	for (n, _, _) in occupied_modes(b)	179✔
410	result += n * (n - 1)	2,054✔
411	end	4,045✔
412	return result	103✔
413	end
414
415	@inline function bose_hubbard_interaction(::Val{1}, b::BoseFS{<:Any,<:Any,<:BitString})	47,151,949✔
416	# currently this ammounts to counting occupation numbers of modes
417	chunk = chunks(b.bs)[1]	47,115,681✔
418	matrixelementint = 0	47,038,381✔
419	while !iszero(chunk)	235,051,101✔
420	chunk >>>= (trailing_zeros(chunk) % UInt) # proceed to next occupied mode	189,230,522✔
421	bosonnumber = trailing_ones(chunk) # count how many bosons inside	189,680,633✔
422	# surpsingly it is faster to not check whether this is nonzero and do the
423	# following operations anyway
424	chunk >>>= (bosonnumber % UInt) # remove the counted mode	189,441,870✔
425	matrixelementint += bosonnumber * (bosonnumber - 1)	189,241,318✔
426	end	189,577,665✔
427	return matrixelementint	46,934,477✔
428	end

RimuQMC / Rimu.jl / 12429834995

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous