Merge branch 'develop' into feature/pdvec

.github/workflows/actions.yml

@@ -58,7 +58,7 @@ jobs:
       - name: "Process coverage"
         uses: julia-actions/julia-processcoverage@latest
       - name: "Coveralls"
-        uses: coverallsapp/github-action@v1.1.2
+        uses: coverallsapp/github-action@v2
         with:
           github-token: ${{ secrets.github_token }}
           path-to-lcov: ./lcov.info

docs/src/hamiltonians.md

@@ -66,6 +66,7 @@ HubbardMom1DEP
 ```@docs
 HOCartesianContactInteractions
 HOCartesianEnergyConservedPerDim
+HOCartesianCentralImpurity
 ```
 
 ### Other
@@ -166,11 +167,12 @@ and [`HOCartesianEnergyConservedPerDim`](@ref).
 ```@docs
 get_all_blocks
 fock_to_cart
-parity_block_seed_addresses
 ```
 Underlying integrals for the interaction matrix elements are implemented in the following unexported functions
 ```@docs
 Hamiltonians.four_oscillator_integral_general
+Hamiltonians.ho_delta_potential
+Hamiltonians.log_abs_oscillator_zero
 ```
 
 ## Index

src/BitStringAddresses/sortedparticlelist.jl

@@ -44,6 +44,15 @@ function SortedParticleList{N,M}() where {N,M}
     T = select_int_type(M)
     return SortedParticleList{N,M,T}(ones(SVector{N,T}))
 end
+function SortedParticleList(onr)
+    N = sum(onr)
+    M = length(onr)
+    T = select_int_type(M)
+    eltype(onr) <: Integer || throw(ArgumentError("the onr can only contain integers"))
+    any(≥(0), onr) || throw(ArgumentError("all elements of the onr should be positive"))
+
+    return from_onr(SortedParticleList{N,M,T}, onr)
+end
 function from_onr(::Type{S}, onr) where {N,M,T,S<:SortedParticleList{N,M,T}}
     spl = zeros(MVector{N,T})
     curr = 1

src/Hamiltonians/HOCartesianCentralImpurity.jl

@@ -0,0 +1,210 @@
+"""
+    log_abs_oscillator_zero(n)
+
+Compute the logarithm of the absolute value of the ``n^\\mathrm{th}`` 1D 
+harmonic oscillator function evaluated at the origin. The overall sign is
+determined when the matrix element is evaluated in [`ho_delta_potential`](@ref).
+"""
+function log_abs_oscillator_zero(n)
+    isodd(n) && return -Inf, 1.0
+    x, _ = SpecialFunctions.logabsgamma((1-n)/2)
+    y, _ = SpecialFunctions.logabsgamma(n+1)
+    result = log(pi)/4 + (n/2)*log(2) - x - y/2
+    return result
+end
+
+"""
+    ho_delta_potential(S, i, j; [vals])
+
+Returns the matrix element of a delta potential at the centre of a trap, i.e.
+the  product of two harmonic oscillator functions evaluated at the origin, 
+```math
+    v_{ij} = \\phi_{\\mathbf{n}_i}(0) \\phi_{\\mathbf{n}_j}(0)
+```
+which is only non-zero for even-parity states. The `i`th single particle state 
+corresponds to a ``D``-tuple of harmonic oscillator indices ``\\mathbf{n}_i``. 
+`S` defines the bounds of Cartesian harmonic oscillator indices for each dimension.
+The optional keyword argument `vals` allows passing pre-computed values of 
+``\\phi_i(0)`` to speed-up the calculation. The values can be calculated with 
+[`log_abs_oscillator_zero`](@ref).
+
+See also [`HOCartesianCentralImpurity`](@ref).
+"""
+function ho_delta_potential(S, i, j; 
+    vals = [log_abs_oscillator_zero(k) for k in 0:2:(maximum(S)-1)]
+    )
+    states = CartesianIndices(S)
+    ho_indices = (Tuple(states[i])..., Tuple(states[j])...) .- 1
+    if all(iseven.(ho_indices))
+        doubly_evens = count(iseven.(ho_indices .÷ 2))
+        sign = (-1)^doubly_evens
+        return sign * exp(sum(vals[k÷2 + 1] for k in ho_indices))
+    else
+        return 0.0
+    end
+end
+
+"""
+    HOCartesianCentralImpurity(addr; kwargs...)
+
+Hamiltonian of non-interacting particles in an arbitrary harmonic trap with a delta-function 
+potential at the centre, with strength `g`,
+```math
+\\hat{H}_\\mathrm{rel} = \\sum_\\mathbf{i} ϵ_\\mathbf{i} n_\\mathbf{i} 
+    + g\\sum_\\mathbf{ij} V_\\mathbf{ij} a^†_\\mathbf{i} a_\\mathbf{j}.
+```
+For a ``D``-dimensional harmonic oscillator indices ``\\mathbf{i}, \\mathbf{j}, \\ldots`` 
+are ``D``-tuples. The energy scale is defined by the first dimension i.e. ``\\hbar \\omega_x`` 
+so that single particle energies are 
+```math
+    \\frac{\\epsilon_\\mathbf{i}}{\\hbar \\omega_x} = (i_x + 1/2) + \\eta_y (i_y+1/2) + \\ldots.
+```
+The factors ``\\eta_y, \\ldots`` allow for anisotropic trapping geometries and are assumed to 
+be greater than `1` so that ``\\omega_x`` is the smallest trapping frequency.
+
+Matrix elements ``V_{\\mathbf{ij}}`` are for a delta function potential calculated in this basis
+```math
+    V_{\\mathbf{ij}} = \\prod_{d \\in x, y,\\ldots} \\psi_{i_d}(0) \\psi_{j_d}(0).
+```
+Only even parity states feel this impurity, so all ``i_d`` are even. Note that the matrix 
+representation of this Hamiltonian for a single particle is completely dense in the even-parity 
+subspace.
+
+# Arguments
+
+* `addr`: the starting address, defines number of particles and total number of modes.
+* `max_nx = num_modes(addr) - 1`: the maximum harmonic oscillator index number in the ``x``-dimension.
+    Must be even. Index number for the harmonic oscillator groundstate is `0`.
+* `ηs = ()`: a tuple of aspect ratios for the remaining dimensions `(η_y, ...)`. Should be empty 
+    for a 1D trap or contain values greater than `1.0`. The maximum index 
+    in other dimensions will be the largest even number less than `M/η_y`.
+* `S = nothing`: Instead of `max_nx`, manually set the number of levels in each dimension, 
+    including the groundstate. Must be a `Tuple` of `Int`s.
+* `g = 1.0`: the strength of the delta impurity in (``x``-dimension) trap units.
+* `impurity_only=false`: if set to `true` then the trap energy terms are ignored. Useful if 
+    only energy shifts due to the impurity are required.
+
+!!! warning
+        Due to use of `SpecialFunctions` with large arguments the matrix representation of 
+        this Hamiltonian may not be strictly symmetric, but is approximately symmetric within
+        machine precision.
+
+See also [`HOCartesianContactInteractions`](@ref) and[`HOCartesianEnergyConservedPerDim`](@ref).
+"""
+struct HOCartesianCentralImpurity{
+    D,  # number of dimensions
+    A
+} <: AbstractHamiltonian{Float64}
+    addr::A
+    S::NTuple{D,Int64}
+    aspect::NTuple{D,Float64}
+    vtable::Vector{Float64}  # interaction coefficients
+    g::Float64
+    impurity_only::Bool
+end
+
+function HOCartesianCentralImpurity(
+        addr::SingleComponentFockAddress; 
+        max_nx::Int = num_modes(addr) - 1,
+        S = nothing,
+        ηs = (), 
+        g = 1.0,
+        impurity_only = false
+    )
+    any(ηs .< 1.0) && throw(ArgumentError("Aspect ratios must be greater than 1.0"))
+    if isnothing(S)
+        iseven(max_nx) || throw(ArgumentError("max_nx must be even"))
+        D = length(ηs) + 1
+        Srest = map(x -> 2floor(Int, max_nx÷2 / x) + 1, ηs)
+        S = (max_nx + 1, Srest...)
+    else
+        D = length(S)
+    end
+    aspect = (1.0, float.(ηs)...)
+
+    num_modes(addr) == prod(S) || throw(ArgumentError("number of modes does not match starting address"))
+
+    v_vec = [log_abs_oscillator_zero(i) for i in 0:2:S[1]-1]
+
+    return HOCartesianCentralImpurity{D,typeof(addr)}(addr, S, aspect, v_vec, g, impurity_only)
+end
+
+function Base.show(io::IO, H::HOCartesianCentralImpurity)
+    print(io, "HOCartesianCentralImpurity($(H.addr); max_nx=$(H.S[1]-1), ηs=$(H.aspect[2:end]), g=$(H.g), impurity_only=$(H.impurity_only))")
+end
+
+Base.:(==)(H::HOCartesianCentralImpurity, G::HOCartesianCentralImpurity) = all(map(p -> getproperty(H, p) == getproperty(G, p), propertynames(H)))
+
+starting_address(H::HOCartesianCentralImpurity) = H.addr
+
+LOStructure(::Type{<:HOCartesianCentralImpurity}) = IsHermitian()
+
+### DIAGONAL ELEMENTS ###
+@inline function noninteracting_energy(S, aspect, omm::BoseOccupiedModeMap)
+    states = CartesianIndices(S)
+    return sum(omm) do p
+        indices = Tuple(states[p.mode]) .- 1
+        dot(indices, aspect) + sum(aspect)/2
+    end
+end
+noninteracting_energy(H::HOCartesianCentralImpurity, addr) = noninteracting_energy(H.S, H.aspect, OccupiedModeMap(addr))
+
+@inline function diagonal_element(H::HOCartesianCentralImpurity, addr)
+    omm = OccupiedModeMap(addr)
+    u = H.g / sqrt(prod(H.aspect))
+    result = u * sum(p -> ho_delta_potential(H.S, p.mode, p.mode; vals = H.vtable), omm)
+    if !H.impurity_only
+        result += noninteracting_energy(H, addr)
+    end
+    return result
+end
+
+### OFFDIAGONAL ELEMENTS ###
+get_offdiagonal(H::HOCartesianCentralImpurity, addr, i) = offdiagonals(H, addr)[i]
+num_offdiagonals(H::HOCartesianCentralImpurity, addr) = (num_modes(addr) - 1) * length(occupied_modes(addr))
+
+###
+### offdiagonals
+###
+"""
+    HOCartImpurityOffdiagonals
+
+Specialized [`AbstractOffdiagonals`](@ref) for [`HOCartesianCentralImpurity`](@ref).
+"""
+struct HOCartImpurityOffdiagonals{
+    A<:BoseFS,O<:OccupiedModeMap,H<:HOCartesianCentralImpurity
+} <: AbstractOffdiagonals{A,Float64}
+    ham::H
+    address::A
+    omm::O
+    u::Float64
+    length::Int    
+end
+
+function offdiagonals(H::HOCartesianCentralImpurity, addr::SingleComponentFockAddress)
+    omm = OccupiedModeMap(addr)
+    num_offs = num_offdiagonals(H, addr)
+    u = H.g / sqrt(prod(H.aspect))
+    return HOCartImpurityOffdiagonals(H, addr, omm, u, num_offs)
+end
+
+function Base.getindex(offs::HOCartImpurityOffdiagonals, chosen)
+    addr = offs.address
+    omm = offs.omm
+    S = offs.ham.S
+    vals = offs.ham.vtable    
+    u = offs.u
+
+    P = num_modes(addr)
+    i, j = fldmod1(chosen, P - 1)
+    index_i = omm[i]
+    j += j ≥ index_i.mode   # skip self move
+    index_j = find_mode(addr, j)
+    new_addr, val = excitation(addr, (index_j,), (index_i,))
+
+    impurity = ho_delta_potential(S, index_i.mode, index_j.mode; vals)
+
+    return new_addr, val * impurity * u
+end
+
+Base.size(s::HOCartImpurityOffdiagonals) = (s.length,)

src/Hamiltonians/HOCartesianEnergyConservedPerDim.jl

@@ -302,8 +302,7 @@ end
 """
     HOCartSeparableOffdiagonals
 
-Specialized [`AbstractOffdiagonals`](@ref) that keeps track of singly and doubly occupied
-sites in current address.
+Specialized [`AbstractOffdiagonals`](@ref) for [`HOCartesianEnergyConservedPerDim`](@ref).
 """
 struct HOCartSeparableOffdiagonals{
     A<:BoseFS,T,H<:AbstractHamiltonian{T},O<:OccupiedModeMap

src/Hamiltonians/Hamiltonians.jl

@@ -21,6 +21,7 @@ Momentum space Hubbard models
 Harmonic oscillator models
 - [`HOCartesianContactInteractions`](@ref)
 - [`HOCartesianEnergyConservedPerDim`](@ref)
+- [`HOCartesianCentralImpurity`](@ref)
 
 Other
 - [`MatrixHamiltonian`](@ref)
@@ -80,8 +81,9 @@ export num_neighbours, neighbour_site, num_dimensions
 
 export sparse # from SparseArrays
 
-export HOCartesianContactInteractions, HOCartesianEnergyConservedPerDim, AxialAngularMomentumHO
-export get_all_blocks, fock_to_cart, parity_block_seed_addresses
+export HOCartesianContactInteractions, HOCartesianEnergyConservedPerDim, HOCartesianCentralImpurity 
+export AxialAngularMomentumHO
+export get_all_blocks, fock_to_cart
 
 include("abstract.jl")
 include("offdiagonals.jl")
@@ -114,6 +116,7 @@ include("Momentum.jl")
 
 include("HOCartesianContactInteractions.jl")
 include("HOCartesianEnergyConservedPerDim.jl")
+include("HOCartesianCentralImpurity.jl")
 include("vertices.jl")
 include("ho-cart-tools.jl")
 include("angular_momentum.jl")

src/Hamiltonians/ho-cart-tools.jl

@@ -29,10 +29,9 @@ Keyword arguments:
     Useful for anisotropic traps
 
 Note: If `h` was constructed with option `block_by_level = false` then the block seeds 
-`addr` are determined by parity, see [`parity_block_seed_addresses`](@ref). In this case 
-the blocks are not generated; `t_basis` will be zero, and `block_size` will be an 
-estimate. Pass the seed addresses to [`BasisSetRep`](@ref) with an appropriate `filter` 
-to generate the blocks.
+`addr` are determined by parity. In this case the blocks are not generated; `t_basis` 
+will be zero, and `block_size` will be an estimate. Pass the seed addresses to 
+[`BasisSetRep`](@ref) with an appropriate `filter` to generate the blocks.
 """
 function get_all_blocks(h::Union{HOCartesianContactInteractions{D,<:Any,B},HOCartesianEnergyConservedPerDim{D}}; 
         target_energy = nothing, 
@@ -171,6 +170,29 @@ function get_all_blocks_comb(h;
     return df
 end
 
+"""
+    parity_block_seed_addresses(h::HOCartesianContactInteractions{D})
+
+Get a vector of addresses that each have different parity with respect to 
+the trap geometry defined by the Hamiltonian `H`. The result will have `2^D`
+`BoseFS` addresses for a `D`-dimensional trap. This is useful for 
+[`HOCartesianContactInteractions`](@ref) with option `block_by_level = false`.
+"""
+function parity_block_seed_addresses(h::HOCartesianContactInteractions{D,A}) where {D,A}
+    P = prod(h.S)
+    N = num_particles(h.addr)
+    breakpoints = accumulate(*, (1, h.S[1:end-1]...))
+    seeds = A[]
+    for c in with_replacement_combinations([0,1], D), p in multiset_permutations(c, D)
+        index = 1 + dot(p, breakpoints)
+        push!(seeds,
+            BoseFS(P, index => 1, 1 => N-1)
+        )
+    end
+
+    return seeds
+end
+
 # specialised version if `block_by_level = false`
 function get_all_blocks_parity(h::HOCartesianContactInteractions{D,<:Any,B}) where {D,B}
     # check if H blocks by level
@@ -207,27 +229,4 @@ function fock_to_cart(
         OccupiedModeMap(addr))...)
 
     return SVector{N,NTuple{D,Int}}(cart)
-end
-
-"""
-    parity_block_seed_addresses(h::HOCartesianContactInteractions{D})
-
-Get a vector of addresses that each have different parity with respect to 
-the trap geometry defined by the Hamiltonian `H`. The result will have `2^D`
-`BoseFS` addresses for a `D`-dimensional trap. This is useful for 
-[`HOCartesianContactInteractions`](@ref) with option `block_by_level = false`.
-"""
-function parity_block_seed_addresses(h::HOCartesianContactInteractions{D,A}) where {D,A}
-    P = prod(h.S)
-    N = num_particles(h.addr)
-    breakpoints = accumulate(*, (1, h.S[1:end-1]...))
-    seeds = A[]
-    for c in with_replacement_combinations([0,1], D), p in multiset_permutations(c, D)
-        index = 1 + dot(p, breakpoints)
-        push!(seeds,
-            BoseFS(P, index => 1, 1 => N-1)
-        )
-    end
-
-    return seeds
 end

src/Rimu.jl

@@ -43,7 +43,7 @@ include("StochasticStyles/StochasticStyles.jl")
 @reexport using .StochasticStyles
 include("DictVectors/DictVectors.jl")
 @reexport using .DictVectors
-include("RimuIO.jl")
+include("RimuIO/RimuIO.jl")
 @reexport using .RimuIO
 include("StatsTools/StatsTools.jl")
 @reexport using .StatsTools