Skip to content

Commit

Permalink
compact printing of addresses in Hamiltonians
Browse files Browse the repository at this point in the history
  • Loading branch information
joachimbrand committed May 24, 2024
1 parent a649be5 commit 049a634
Show file tree
Hide file tree
Showing 15 changed files with 191 additions and 185 deletions.
9 changes: 4 additions & 5 deletions src/ExactDiagonalization/basis_set_representation.jl
Original file line number Diff line number Diff line change
Expand Up @@ -188,17 +188,16 @@ are passed on to `Base.sortperm`.
julia> hamiltonian = HubbardReal1D(BoseFS(1,0,0));
julia> bsr = BasisSetRepresentation(hamiltonian)
BasisSetRepresentation(HubbardReal1D(BoseFS{1,3}(1, 0, 0); u=1.0, t=1.0)) with dimension 3 and 9 stored entries:3×3 SparseArrays.SparseMatrixCSC{Float64, Int64} with 9 stored entries:
BasisSetRepresentation(HubbardReal1D(fs"|1 0 0⟩"; u=1.0, t=1.0)) with dimension 3 and 9 stored entries:3×3 SparseArrays.SparseMatrixCSC{Float64, Int64} with 9 stored entries:
0.0 -1.0 -1.0
-1.0 0.0 -1.0
-1.0 -1.0 0.0
julia> BasisSetRepresentation(hamiltonian, bsr.basis[1:2]; filter = Returns(false)) # passing addresses and truncating
BasisSetRepresentation(HubbardReal1D(BoseFS{1,3}(1, 0, 0); u=1.0, t=1.0)) with dimension 2 and 4 stored entries:2×2 SparseArrays.SparseMatrixCSC{Float64, Int64} with 4 stored entries:
BasisSetRepresentation(HubbardReal1D(fs"|1 0 0⟩"; u=1.0, t=1.0)) with dimension 2 and 4 stored entries:2×2 SparseArrays.SparseMatrixCSC{Float64, Int64} with 4 stored entries:
0.0 -1.0
-1.0 0.0
```
```julia-repl
julia> using LinearAlgebra; eigvals(Matrix(bsr)) # eigenvalues
3-element Vector{Float64}:
-1.9999999999999996
Expand All @@ -216,7 +215,7 @@ DVec{BoseFS{1, 3, BitString{3, 1, UInt8}},Float64} with 3 entries, style = IsDet
fs"|0 0 1⟩" => -0.57735
fs"|0 1 0⟩" => -0.57735
fs"|1 0 0⟩" => -0.57735
```
```
Has methods for [`dimension`](@ref), [`sparse`](@ref), [`Matrix`](@ref),
[`starting_address`](@ref).
Expand Down
36 changes: 18 additions & 18 deletions src/Hamiltonians/BoseHubbardMom1D2C.jl
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
"""
BoseHubbardMom1D2C(add::BoseFS2C; ua=1.0, ub=1.0, ta=1.0, tb=1.0, v=1.0, kwargs...)
BoseHubbardMom1D2C(address::BoseFS2C; ua=1.0, ub=1.0, ta=1.0, tb=1.0, v=1.0, kwargs...)
Implements a one-dimensional Bose Hubbard chain in momentum space with a two-component
Bose gas.
Expand All @@ -10,7 +10,7 @@ Bose gas.
# Arguments
* `add`: the starting address.
* `address`: the starting address.
* `ua`: the `u` parameter for Hamiltonian a.
* `ub`: the `u` parameter for Hamiltonian b.
* `ta`: the `t` parameter for Hamiltonian a.
Expand All @@ -29,9 +29,9 @@ struct BoseHubbardMom1D2C{T,HA,HB,V} <: TwoComponentHamiltonian{T}
hb::HB
end

function BoseHubbardMom1D2C(add::BoseFS2C; ua=1.0, ub=1.0, ta=1.0, tb=1.0, v=1.0, args...)
ha = HubbardMom1D(add.bsa; u=ua, t=ta, args...)
hb = HubbardMom1D(add.bsb; u=ub, t=tb, args...)
function BoseHubbardMom1D2C(address::BoseFS2C; ua=1.0, ub=1.0, ta=1.0, tb=1.0, v=1.0, args...)
ha = HubbardMom1D(address.bsa; u=ua, t=ta, args...)
hb = HubbardMom1D(address.bsb; u=ub, t=tb, args...)
T = promote_type(eltype(ha), eltype(hb))
return BoseHubbardMom1D2C{T,typeof(ha),typeof(hb),v}(ha, hb)
end
Expand Down Expand Up @@ -59,18 +59,18 @@ Base.getproperty(h::BoseHubbardMom1D2C, ::Val{:ha}) = getfield(h, :ha)
Base.getproperty(h::BoseHubbardMom1D2C, ::Val{:hb}) = getfield(h, :hb)
Base.getproperty(h::BoseHubbardMom1D2C{<:Any,<:Any,<:Any,V}, ::Val{:v}) where {V} = V

function num_offdiagonals(ham::BoseHubbardMom1D2C, add::BoseFS2C)
M = num_modes(add)
sa = num_occupied_modes(add.bsa)
sb = num_occupied_modes(add.bsb)
return num_offdiagonals(ham.ha, add.bsa) + num_offdiagonals(ham.hb, add.bsb) + sa*(M-1)*sb
function num_offdiagonals(ham::BoseHubbardMom1D2C, address::BoseFS2C)
M = num_modes(address)
sa = num_occupied_modes(address.bsa)
sb = num_occupied_modes(address.bsb)
return num_offdiagonals(ham.ha, address.bsa) + num_offdiagonals(ham.hb, address.bsb) + sa*(M-1)*sb
# number of excitations that can be made
end

function diagonal_element(ham::BoseHubbardMom1D2C, add::BoseFS2C)
M = num_modes(add)
onrep_a = onr(add.bsa)
onrep_b = onr(add.bsb)
function diagonal_element(ham::BoseHubbardMom1D2C, address::BoseFS2C)
M = num_modes(address)
onrep_a = onr(address.bsa)
onrep_b = onr(address.bsb)
interaction2c = Int32(0)
for p in 1:M
iszero(onrep_b[p]) && continue
Expand All @@ -79,14 +79,14 @@ function diagonal_element(ham::BoseHubbardMom1D2C, add::BoseFS2C)
end
end
return (
diagonal_element(ham.ha, add.bsa) +
diagonal_element(ham.hb, add.bsb) +
diagonal_element(ham.ha, address.bsa) +
diagonal_element(ham.hb, address.bsb) +
ham.v/M*interaction2c
)
end

function get_offdiagonal(ham::BoseHubbardMom1D2C, add::BoseFS2C, chosen)
return offdiagonals(ham, add)[chosen]
function get_offdiagonal(ham::BoseHubbardMom1D2C, address::BoseFS2C, chosen)
return offdiagonals(ham, address)[chosen]
end

"""
Expand Down
30 changes: 15 additions & 15 deletions src/Hamiltonians/BoseHubbardReal1D2C.jl
Original file line number Diff line number Diff line change
Expand Up @@ -27,9 +27,9 @@ struct BoseHubbardReal1D2C{T,HA,HB,V} <: TwoComponentHamiltonian{T}
hb::HB
end

function BoseHubbardReal1D2C(add::BoseFS2C; ua=1.0,ub=1.0,ta=1.0,tb=1.0,v=1.0)
ha = HubbardReal1D(add.bsa; u=ua, t=ta)
hb = HubbardReal1D(add.bsb; u=ub, t=tb)
function BoseHubbardReal1D2C(address::BoseFS2C; ua=1.0,ub=1.0,ta=1.0,tb=1.0,v=1.0)
ha = HubbardReal1D(address.bsa; u=ua, t=ta)
hb = HubbardReal1D(address.bsb; u=ub, t=tb)
T = promote_type(eltype(ha), eltype(hb))
return BoseHubbardReal1D2C{T,typeof(ha),typeof(hb),v}(ha, hb)
end
Expand Down Expand Up @@ -58,18 +58,18 @@ Base.getproperty(h::BoseHubbardReal1D2C, ::Val{:hb}) = getfield(h, :hb)
Base.getproperty(h::BoseHubbardReal1D2C{<:Any,<:Any,<:Any,V}, ::Val{:v}) where {V} = V

# number of excitations that can be made
function num_offdiagonals(ham::BoseHubbardReal1D2C, add)
return 2*(num_occupied_modes(add.bsa) + num_occupied_modes(add.bsb))
function num_offdiagonals(::BoseHubbardReal1D2C, address)
return 2*(num_occupied_modes(address.bsa) + num_occupied_modes(address.bsb))
end

"""
bose_hubbard_2c_interaction(::BoseFS2C)
Compute the interaction between the two components.
"""
function bose_hubbard_2c_interaction(add::BoseFS2C)
c1 = onr(add.bsa)
c2 = onr(add.bsb)
function bose_hubbard_2c_interaction(address::BoseFS2C)
c1 = onr(address.bsa)
c2 = onr(address.bsb)
interaction = zero(eltype(c1))
for site = 1:length(c1)
if !iszero(c2[site])
Expand All @@ -87,18 +87,18 @@ function diagonal_element(ham::BoseHubbardReal1D2C, address::BoseFS2C)
)
end

function get_offdiagonal(ham::BoseHubbardReal1D2C, add, chosen)
nhops = num_offdiagonals(ham,add)
nhops_a = 2 * num_occupied_modes(add.bsa)
function get_offdiagonal(ham::BoseHubbardReal1D2C, address, chosen)
nhops = num_offdiagonals(ham,address)
nhops_a = 2 * num_occupied_modes(address.bsa)
if chosen nhops_a
naddress_from_bsa, onproduct = hopnextneighbour(add.bsa, chosen)
naddress_from_bsa, onproduct = hopnextneighbour(address.bsa, chosen)
elem = - ham.ha.t * onproduct
return BoseFS2C(naddress_from_bsa,add.bsb), elem
return BoseFS2C(naddress_from_bsa,address.bsb), elem
else
chosen -= nhops_a
naddress_from_bsb, onproduct = hopnextneighbour(add.bsb, chosen)
naddress_from_bsb, onproduct = hopnextneighbour(address.bsb, chosen)
elem = -ham.hb.t * onproduct
return BoseFS2C(add.bsa,naddress_from_bsb), elem
return BoseFS2C(address.bsa,naddress_from_bsb), elem
end
# return new address and matrix element
end
Expand Down
13 changes: 7 additions & 6 deletions src/Hamiltonians/ExtendedHubbardReal1D.jl
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,7 @@ Implements the extended Hubbard model on a one-dimensional chain in real space.
"""
struct ExtendedHubbardReal1D{TT,A<:SingleComponentFockAddress,U,V,T} <: AbstractHamiltonian{TT}
add::A
address::A
end

# addr for compatibility.
Expand All @@ -26,19 +26,20 @@ function ExtendedHubbardReal1D(addr; u=1.0, v=1.0, t=1.0)
end

function Base.show(io::IO, h::ExtendedHubbardReal1D)
print(io, "ExtendedHubbardReal1D($(h.add); u=$(h.u), v=$(h.v), t=$(h.t))")
compact_addr = repr(h.address, context=:compact => true) # compact print address
print(io, "ExtendedHubbardReal1D($(compact_addr); u=$(h.u), v=$(h.v), t=$(h.t))")
end

function starting_address(h::ExtendedHubbardReal1D)
return getfield(h, :add)
return getfield(h, :address)
end

dimension(::ExtendedHubbardReal1D, address) = number_conserving_dimension(address)

LOStructure(::Type{<:ExtendedHubbardReal1D{<:Real}}) = IsHermitian()

Base.getproperty(h::ExtendedHubbardReal1D, s::Symbol) = getproperty(h, Val(s))
Base.getproperty(h::ExtendedHubbardReal1D, ::Val{:add}) = getfield(h, :add)
Base.getproperty(h::ExtendedHubbardReal1D, ::Val{:address}) = getfield(h, :address)
Base.getproperty(h::ExtendedHubbardReal1D{<:Any,<:Any,U}, ::Val{:u}) where U = U
Base.getproperty(h::ExtendedHubbardReal1D{<:Any,<:Any,<:Any,V}, ::Val{:v}) where V = V
Base.getproperty(h::ExtendedHubbardReal1D{<:Any,<:Any,<:Any,<:Any,T}, ::Val{:t}) where T = T
Expand Down Expand Up @@ -78,7 +79,7 @@ function diagonal_element(h::ExtendedHubbardReal1D, b::SingleComponentFockAddres
return h.u * bhinteraction / 2 + h.v * ebhinteraction
end

function get_offdiagonal(h::ExtendedHubbardReal1D, add::SingleComponentFockAddress, chosen)
naddress, onproduct = hopnextneighbour(add, chosen)
function get_offdiagonal(h::ExtendedHubbardReal1D, address::SingleComponentFockAddress, chosen)
naddress, onproduct = hopnextneighbour(address, chosen)
return naddress, - h.t * onproduct
end
4 changes: 2 additions & 2 deletions src/Hamiltonians/GutzwillerSampling.jl
Original file line number Diff line number Diff line change
Expand Up @@ -21,10 +21,10 @@ After construction, we can access the underlying Hamiltonian with `G.hamiltonian
```jldoctest
julia> H = HubbardMom1D(BoseFS(1,1,1); u=6.0, t=1.0)
HubbardMom1D(BoseFS{3,3}(1, 1, 1); u=6.0, t=1.0)
HubbardMom1D(fs"|1 1 1⟩"; u=6.0, t=1.0)
julia> G = GutzwillerSampling(H, g=0.3)
GutzwillerSampling(HubbardMom1D(BoseFS{3,3}(1, 1, 1); u=6.0, t=1.0); g=0.3)
GutzwillerSampling(HubbardMom1D(fs"|1 1 1⟩"; u=6.0, t=1.0); g=0.3)
julia> get_offdiagonal(H, BoseFS(2, 1, 0), 1)
(BoseFS{3,3}(1, 0, 2), 2.0)
Expand Down
53 changes: 27 additions & 26 deletions src/Hamiltonians/HOCartesianCentralImpurity.jl
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
"""
log_abs_oscillator_zero(n)
Compute the logarithm of the absolute value of the ``n^\\mathrm{th}`` 1D
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).
"""
Expand All @@ -17,20 +17,20 @@ 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,
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``.
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
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;
function ho_delta_potential(S, i, j;
vals = [log_abs_oscillator_zero(k) for k in 0:2:(maximum(S)-1)]
)
states = CartesianIndices(S)
Expand All @@ -47,45 +47,45 @@ end
"""
HOCartesianCentralImpurity(addr; kwargs...)
Hamiltonian of non-interacting particles in an arbitrary harmonic trap with a delta-function
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}
\\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
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
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
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
* `η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,
* `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
* `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
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.
Expand All @@ -104,10 +104,10 @@ struct HOCartesianCentralImpurity{
end

function HOCartesianCentralImpurity(
addr::SingleComponentFockAddress;
addr::SingleComponentFockAddress;
max_nx::Int = num_modes(addr) - 1,
S = nothing,
ηs = (),
ηs = (),
g = 1.0,
impurity_only = false
)
Expand All @@ -121,7 +121,7 @@ function HOCartesianCentralImpurity(
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]
Expand All @@ -130,7 +130,8 @@ function HOCartesianCentralImpurity(
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))")
compact_addr = repr(H.addr, context=:compact => true) # compact print address
print(io, "HOCartesianCentralImpurity($(compact_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)))
Expand Down Expand Up @@ -178,7 +179,7 @@ struct HOCartImpurityOffdiagonals{
address::A
omm::O
u::Float64
length::Int
length::Int
end

function offdiagonals(H::HOCartesianCentralImpurity, addr::SingleComponentFockAddress)
Expand All @@ -192,7 +193,7 @@ function Base.getindex(offs::HOCartImpurityOffdiagonals, chosen)
addr = offs.address
omm = offs.omm
S = offs.ham.S
vals = offs.ham.vtable
vals = offs.ham.vtable
u = offs.u

P = num_modes(addr)
Expand Down
Loading

0 comments on commit 049a634

Please sign in to comment.