diff --git a/src/ExactDiagonalization/basis_set_representation.jl b/src/ExactDiagonalization/basis_set_representation.jl
index 9b7669804..f9f8bf81f 100644
--- a/src/ExactDiagonalization/basis_set_representation.jl
+++ b/src/ExactDiagonalization/basis_set_representation.jl
@@ -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
@@ -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).
 
diff --git a/src/Hamiltonians/BoseHubbardMom1D2C.jl b/src/Hamiltonians/BoseHubbardMom1D2C.jl
index 181b783db..76c12c91e 100644
--- a/src/Hamiltonians/BoseHubbardMom1D2C.jl
+++ b/src/Hamiltonians/BoseHubbardMom1D2C.jl
@@ -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.
@@ -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.
@@ -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
@@ -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
@@ -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
 
 """
diff --git a/src/Hamiltonians/BoseHubbardReal1D2C.jl b/src/Hamiltonians/BoseHubbardReal1D2C.jl
index 9a4aaebe8..74411e6d4 100644
--- a/src/Hamiltonians/BoseHubbardReal1D2C.jl
+++ b/src/Hamiltonians/BoseHubbardReal1D2C.jl
@@ -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
@@ -58,8 +58,8 @@ 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
 
 """
@@ -67,9 +67,9 @@ end
 
 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])
@@ -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
diff --git a/src/Hamiltonians/ExtendedHubbardReal1D.jl b/src/Hamiltonians/ExtendedHubbardReal1D.jl
index 47bb3d6ee..7079497d7 100644
--- a/src/Hamiltonians/ExtendedHubbardReal1D.jl
+++ b/src/Hamiltonians/ExtendedHubbardReal1D.jl
@@ -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.
@@ -26,11 +26,12 @@ 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)
@@ -38,7 +39,7 @@ dimension(::ExtendedHubbardReal1D, address) = number_conserving_dimension(addres
 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
@@ -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
diff --git a/src/Hamiltonians/GutzwillerSampling.jl b/src/Hamiltonians/GutzwillerSampling.jl
index b5c88cf80..fd3f7bd66 100644
--- a/src/Hamiltonians/GutzwillerSampling.jl
+++ b/src/Hamiltonians/GutzwillerSampling.jl
@@ -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)
diff --git a/src/Hamiltonians/HOCartesianCentralImpurity.jl b/src/Hamiltonians/HOCartesianCentralImpurity.jl
index 6223d0251..0e48ff46d 100644
--- a/src/Hamiltonians/HOCartesianCentralImpurity.jl
+++ b/src/Hamiltonians/HOCartesianCentralImpurity.jl
@@ -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).
 """
@@ -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)
@@ -47,27 +47,27 @@ 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
@@ -75,17 +75,17 @@ subspace.
 * `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.
 
@@ -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
     )
@@ -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]
@@ -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)))
@@ -178,7 +179,7 @@ struct HOCartImpurityOffdiagonals{
     address::A
     omm::O
     u::Float64
-    length::Int    
+    length::Int
 end
 
 function offdiagonals(H::HOCartesianCentralImpurity, addr::SingleComponentFockAddress)
@@ -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)
diff --git a/src/Hamiltonians/HOCartesianContactInteractions.jl b/src/Hamiltonians/HOCartesianContactInteractions.jl
index bb7d3d47e..6bcc77926 100644
--- a/src/Hamiltonians/HOCartesianContactInteractions.jl
+++ b/src/Hamiltonians/HOCartesianContactInteractions.jl
@@ -1,18 +1,18 @@
 """
     four_oscillator_integral_general(i, j, k, l; max_level = typemax(Int))
 
-Integral of four one-dimensional harmonic oscillator functions, 
+Integral of four one-dimensional harmonic oscillator functions,
 ```math
-    \\mathcal{I}(i,j,k,l) = \\int_{-\\infty}^\\infty dx \\, 
+    \\mathcal{I}(i,j,k,l) = \\int_{-\\infty}^\\infty dx \\,
     \\phi_i(x) \\phi_j(x) \\phi_k(x) \\phi_l(x)
 ```
 Indices `i,j,k,l` start at `0` for the groundstate.
 
-This integral has a closed form in terms of the hypergeometric ``_{3}F_2`` function, 
-and is non-zero unless ``i+j+k+l`` is odd. See e.g. 
+This integral has a closed form in terms of the hypergeometric ``_{3}F_2`` function,
+and is non-zero unless ``i+j+k+l`` is odd. See e.g.
 [Titchmarsh (1948)](https://doi.org/10.1112/jlms/s1-23.1.15).
-This is a generalisation of the closed form in 
-[Papenbrock (2002)](https://doi.org/10.1103/PhysRevA.65.033606), which is is the special 
+This is a generalisation of the closed form in
+[Papenbrock (2002)](https://doi.org/10.1103/PhysRevA.65.033606), which is is the special
 case where ``i+j == k+l``, but is numerically unstable for large arguments.
 Used in [`HOCartesianContactInteractions`](@ref) and [`HOCartesianEnergyConservedPerDim`](@ref).
 """
@@ -22,15 +22,15 @@ function four_oscillator_integral_general(i, j, k, l; max_level = typemax(Int))
 
     # enforce correct symmetry and k ≥ l
     i, j, k, l = sort(SVector(i, j, k, l); rev = true)
-    
-    a = i + j - k + l + 1    
+
+    a = i + j - k + l + 1
     b = i - j + k - l + 1
     c = -i + j + k - l + 1
     d = -i - j + k - l + 1
-    
+
     p = sqrt(2 * gamma(i + 1) * gamma(j + 1) * gamma(k + 1) * gamma(l + 1)) * pi^2
     q = gamma(a/2) * gamma(b/2) * gamma(c/2) * gamma(k + 1) / gamma(k - l + 1)
-    
+
     f1 = _₃F₂(float(-l), b/2, c/2, float(1 + k - l), d/2, 1.0)
 
     if isnan(f1)
@@ -41,14 +41,14 @@ function four_oscillator_integral_general(i, j, k, l; max_level = typemax(Int))
     else
         f = f1
     end
-    
+
     return f * q / p
 end
 
 """
     largest_two_point_interval(i, j, [v=1,] w) -> range
 
-For two points `i` and `j` on a range `v:w`, find the largest subinterval of sites 
+For two points `i` and `j` on a range `v:w`, find the largest subinterval of sites
 such that moving `i` to one of those sites and moving `j` by an equal but opposite
 amount leaves both within `v:w`.
 
@@ -81,7 +81,7 @@ function two_point_interval_bounds(i, j, v, w)
     return left_edge, right_edge
 end
 
-# useful conversion of box to aspect ratios; 
+# useful conversion of box to aspect ratios;
 box_to_aspect(S::NTuple{<:Any,Int}) = (S[1] - 1) .// (S .- 1)
 
 @inline function mode_to_energy(i, S, aspect)
@@ -97,11 +97,11 @@ end
 """
     find_Ebounds(i, j, S)
 
-Find the range of single particle energies that a particle in mode `i` can move to while 
+Find the range of single particle energies that a particle in mode `i` can move to while
 preserving total energy of particles `i` and `j`, and constrained to box `S`.
 """
 function find_Ebounds(i, j, S, aspect)
-    
+
     E_i = mode_to_energy(i, S, aspect)
     E_j = mode_to_energy(j, S, aspect)
 
@@ -112,58 +112,58 @@ end
 """
     HOCartesianContactInteractions(addr; S, η, g = 1.0, interaction_only = false, block_by_level = true)
 
-Implements a bosonic harmonic oscillator in Cartesian basis with contact interactions 
+Implements a bosonic harmonic oscillator in Cartesian basis with contact interactions
 ```math
-\\hat{H} = \\sum_{i} \\epsilon_\\mathbf{i} n_\\mathbf{i} + \\frac{g}{2}\\sum_\\mathbf{ijkl} 
+\\hat{H} = \\sum_{i} \\epsilon_\\mathbf{i} n_\\mathbf{i} + \\frac{g}{2}\\sum_\\mathbf{ijkl}
     V_\\mathbf{ijkl} a^†_\\mathbf{i} a^†_\\mathbf{j} a_\\mathbf{k} a_\\mathbf{l}.
 ```
 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 
+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.
 
-By default the offdiagonal elements due to the interactions are consistent with first-order 
+By default the offdiagonal elements due to the interactions are consistent with first-order
 degenerate perturbation theory
 ```math
     V_{\\mathbf{ijkl}} = \\delta_{\\epsilon_\\mathbf{i} + \\epsilon_\\mathbf{j}}
-        ^{\\epsilon_\\mathbf{k} + \\epsilon_\\mathbf{l}} 
+        ^{\\epsilon_\\mathbf{k} + \\epsilon_\\mathbf{l}}
         \\prod_{d \\in x, y,\\ldots} \\mathcal{I}(i_d,j_d,k_d,l_d),
 ```
 where the ``\\delta`` function indicates that the *total* noninteracting energy is conserved
-meaning all states with the same noninteracting energy are connected by this interaction and 
+meaning all states with the same noninteracting energy are connected by this interaction and
 the Hamiltonian blocks according to noninteracting energy levels.
-Setting `block_by_level = false` will disable this restriction and allow coupling between 
-basis states of any noninteracting energy level, leading to many more offdiagonals and 
+Setting `block_by_level = false` will disable this restriction and allow coupling between
+basis states of any noninteracting energy level, leading to many more offdiagonals and
 fewer but larger blocks (the blocks are still distinguished by parity of basis states).
-Alternatively, see [`HOCartesianEnergyConservedPerDim`](@ref) for a model with the stronger 
+Alternatively, see [`HOCartesianEnergyConservedPerDim`](@ref) for a model with the stronger
 restriction that conserves energy separately per spatial dimension.
-The integral ``\\mathcal{I}(a,b,c,d)`` is of four one dimensional harmonic oscillator 
+The integral ``\\mathcal{I}(a,b,c,d)`` is of four one dimensional harmonic oscillator
 basis functions, implemented in [`four_oscillator_integral_general`](@ref).
 
 # Arguments
 
 * `addr`: the starting address, defines number of particles and total number of modes.
-* `S`: Tuple of the number of levels in each dimension, including the groundstate. The 
-    allowed couplings between states is defined by the aspect ratio of `S .- 1`. Defaults 
-    to a 1D spectrum with number of levels matching modes of `addr`. Will be sorted to make 
+* `S`: Tuple of the number of levels in each dimension, including the groundstate. The
+    allowed couplings between states is defined by the aspect ratio of `S .- 1`. Defaults
+    to a 1D spectrum with number of levels matching modes of `addr`. Will be sorted to make
     the first dimension the largest.
 * `η`: Define a custom aspect ratio for the trapping potential strengths, instead of deriving
-    from `S .- 1`. This will only affect the single particle energy scale and not the 
-    interactions. The values are always scaled relative to the first dimension, which sets 
+    from `S .- 1`. This will only affect the single particle energy scale and not the
+    interactions. The values are always scaled relative to the first dimension, which sets
     the energy scale of the system, ``\\hbar\\omega_x``.
-* `g`: the (isotropic) bare interaction parameter. The value of `g` is assumed 
+* `g`: the (isotropic) bare interaction parameter. The value of `g` is assumed
     to be in trap units.
-* `interaction_only`: if set to `true` then the noninteracting single-particle terms are 
+* `interaction_only`: if set to `true` then the noninteracting single-particle terms are
     ignored. Useful if only energy shifts due to interactions are required.
-* `block_by_level`: if set to false will allow the interactions to couple all states without 
+* `block_by_level`: if set to false will allow the interactions to couple all states without
     comparing their noninteracting energy.
 
 !!! warning
-    `num_offdiagonals` is a bad estimate for this Hamiltonian. Take care when building 
+    `num_offdiagonals` is a bad estimate for this Hamiltonian. Take care when building
     a matrix or using QMC methods. Use [`get_all_blocks`](@ref) first then pass option
     `col_hint = block_size` to [`BasisSetRep`](@ref) to safely build the matrix.
 """
@@ -183,9 +183,9 @@ struct HOCartesianContactInteractions{
 end
 
 function HOCartesianContactInteractions(
-        addr::BoseFS; 
+        addr::BoseFS;
         S = (num_modes(addr),),
-        η = nothing, 
+        η = nothing,
         g = 1.0,
         interaction_only = false,
         block_by_level = true
@@ -194,7 +194,7 @@ function HOCartesianContactInteractions(
     D = length(S)
     S_sort = tuple(sort(collect(S); rev = true)...)
     S == S_sort || @warn("dimensions have been reordered")
-    
+
     P = prod(S)
     P == num_modes(addr) || throw(ArgumentError("number of modes does not match starting address"))
 
@@ -231,10 +231,11 @@ function HOCartesianContactInteractions(
 end
 
 function Base.show(io::IO, h::HOCartesianContactInteractions)
+    compact_addr = repr(h.addr, context=:compact => true) # compact print address
     flag = iszero(h.energies)
     # invert the scaling of u parameter
     g = h.u * 2 * sqrt(prod(h.aspect1))
-    print(io, "HOCartesianContactInteractions($(h.addr); S=$(h.S), η=$(h.aspect1), g=$g, interaction_only=$flag)")
+    print(io, "HOCartesianContactInteractions($(compact_addr); S=$(h.S), η=$(h.aspect1), g=$g, interaction_only=$flag)")
 end
 
 Base.:(==)(H::HOCartesianContactInteractions, G::HOCartesianContactInteractions) = all(map(p -> getproperty(H, p) == getproperty(G, p), propertynames(H)))
@@ -250,7 +251,7 @@ LOStructure(::Type{<:HOCartesianContactInteractions}) = IsHermitian()
 function energy_transfer_diagonal(h::HOCartesianContactInteractions{D}, omm::BoseOccupiedModeMap) where {D}
     result = 0.0
     states = CartesianIndices(h.S)    # 1-indexed
-    
+
     for i in eachindex(omm)
         mode_i, occ_i = omm[i].mode, omm[i].occnum
         idx_i = Tuple(states[mode_i])
@@ -264,7 +265,7 @@ function energy_transfer_diagonal(h::HOCartesianContactInteractions{D}, omm::Bos
             idx_j = Tuple(states[mode_j])
             val = 4 * occ_i * occ_j
             result += prod(h.vtable[idx_i[d],idx_j[d],idx_i[d],idx_j[d]] for d in 1:D) * val
-        end        
+        end
     end
     return result * h.u
 end
@@ -317,7 +318,7 @@ Base.size(::HOCartOffdiagonals) = throw(HOCartOffdiagonalsError)
 Base.length(::HOCartOffdiagonals) = throw(HOCartOffdiagonalsError)
 
 # This is the dumb way to loop through valid states.
-# It should take arguments that define where to begin so that the loop can be restarted later 
+# It should take arguments that define where to begin so that the loop can be restarted later
 # Better way would 'jump ahead' when a mode index goes outside the valid energy range.
 # For that I would need dummy indices and probably while loops
 function loop_over_modes(k_start, l_start, S, aspect, Emin, Emax, Etot)
@@ -360,7 +361,7 @@ function loop_over_pairs(S, aspect, pairs, start, block_by_level)
     mode_l = l_start
     if block_by_level
         Es = find_Ebounds(mode_i, mode_j, S, aspect)
-    else 
+    else
         parity_ij = mod(mode_i + mode_j, 2)
     end
     while true
@@ -377,13 +378,13 @@ function loop_over_pairs(S, aspect, pairs, start, block_by_level)
             return pair_index, (mode_i, mode_j, mode_k, mode_l), true
         end
         pair_index += 1
-        # update looping indices        
+        # update looping indices
         p_i, p_j = pairs[pair_index]
         mode_i = p_i.mode
         mode_j = p_j.mode
         if block_by_level
             Es = find_Ebounds(mode_i, mode_j, S, aspect)
-        else 
+        else
             parity_ij = mod(mode_i + mode_j, 2)
         end
         # start with checking lowest possible states in first dimension (x)
@@ -421,7 +422,7 @@ function Base.iterate(off::HOCartOffdiagonals{<:Any,<:Any,B}, iter_state = (1,1,
         states = CartesianIndices(S)    # 1-indexed
         ho_indices = ntuple(x -> Tuple(states[modes[x]]), 4)
         val *= prod(off.ham.vtable[a,b,c,d] for (a,b,c,d) in zip(ho_indices...))
-        # account for j ≤ i and l ≤ k 
+        # account for j ≤ i and l ≤ k
         val *= (1 + (i ≠ j)) * (1 + (k ≠ l)) * off.ham.u
     end
 
diff --git a/src/Hamiltonians/HOCartesianEnergyConservedPerDim.jl b/src/Hamiltonians/HOCartesianEnergyConservedPerDim.jl
index 4b5431f66..75b561e06 100644
--- a/src/Hamiltonians/HOCartesianEnergyConservedPerDim.jl
+++ b/src/Hamiltonians/HOCartesianEnergyConservedPerDim.jl
@@ -2,7 +2,7 @@
 """
     largest_two_point_box(i, j, dims::Tuple) -> ranges
 
-For a box defined by `dims` containing two points `i` and `j`, find the set of all positions 
+For a box defined by `dims` containing two points `i` and `j`, find the set of all positions
 such that shifting `i` to that position with an opposite shift to `j` leaves both points inbounds.
 
 Arguments `i` and `j` are linear indices, `dims` is a `Tuple` defining the inbounds area.
@@ -14,7 +14,7 @@ function largest_two_point_box(i::Int, j::Int, dims::NTuple{D,Int}) where {D}
     state_j = Tuple(cart[j])
 
     ranges = ntuple(d -> largest_two_point_interval(state_i[d], state_j[d], dims[d]), Val(D))
-    
+
     box_size = prod(length.(ranges))
     return ranges, box_size
 end
@@ -23,7 +23,7 @@ end
     find_chosen_pair_moves(omm::OccupiedModeMap, c, S) -> p_i, p_j, c, box_ranges
 
 Find size of valid moves for chosen pair of indices in `omm`. Returns two valid indices `p_i` and `p_j`
-for the initial pair and a tuple of ranges `box_ranges` defining the subbox of valid moves. 
+for the initial pair and a tuple of ranges `box_ranges` defining the subbox of valid moves.
 The index for the chosen move `c` is updated to be valid for this box.
 Arguments are the `OccupiedModeMap` `omm`, the chosen move index `c`
 and the size of the basis grid `S`.
@@ -36,7 +36,7 @@ function find_chosen_pair_moves(omm::OccupiedModeMap, c, S::Tuple)
             if c - box_size < 1
                 return p_i, p_i, box_ranges, c
             else
-                c -= box_size 
+                c -= box_size
             end
         end
         for j in 1:i-1
@@ -45,7 +45,7 @@ function find_chosen_pair_moves(omm::OccupiedModeMap, c, S::Tuple)
             if c - box_size < 1
                 return p_i, p_j, box_ranges, c
             else
-                c -= box_size 
+                c -= box_size
             end
         end
     end
@@ -55,47 +55,47 @@ end
 """
     HOCartesianEnergyConservedPerDim(addr; S, η, g = 1.0, interaction_only = false)
 
-Implements a bosonic harmonic oscillator in Cartesian basis with contact interactions 
+Implements a bosonic harmonic oscillator in Cartesian basis with contact interactions
 ```math
 \\hat{H} = \\sum_{i} ϵ_i n_i + \\frac{g}{2}\\sum_{ijkl} V_{ijkl} a^†_i a^†_j a_k a_l,
 ```
 with the additional restriction that the interactions only couple states with the same
-energy in each dimension separately. See [`HOCartesianContactInteractions`](@ref) for a model that 
+energy in each dimension separately. See [`HOCartesianContactInteractions`](@ref) for a model that
 conserves total energy.
 
-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{ijkl}}`` are for a contact interaction calculated in this basis using 
+Matrix elements ``V_{\\mathbf{ijkl}}`` are for a contact interaction calculated in this basis using
 first-order degenerate perturbation theory.
 ```math
-    V_{\\mathbf{ijkl}} = \\prod_{d \\in x, y,\\ldots} \\mathcal{I}(i_d,j_d,k_d,l_d) 
+    V_{\\mathbf{ijkl}} = \\prod_{d \\in x, y,\\ldots} \\mathcal{I}(i_d,j_d,k_d,l_d)
         \\delta_{i_d + j_d}^{k_d + l_d},
 ```
 where the ``\\delta``-function indicates that the noninteracting energy is conserved along each
 dimension.
-The integral ``\\mathcal{I}(a,b,c,d)`` is of four one dimensional harmonic oscillator 
-basis functions, see [`four_oscillator_integral_general`](@ref), with the additional restriction 
+The integral ``\\mathcal{I}(a,b,c,d)`` is of four one dimensional harmonic oscillator
+basis functions, see [`four_oscillator_integral_general`](@ref), with the additional restriction
 that energy is conserved in each dimension.
 
 # Arguments
 
 * `addr`: the starting address, defines number of particles and total number of modes.
-* `S`: Tuple of the number of levels in each dimension, including the groundstate. Defaults 
-    to a 1D spectrum with number of levels matching modes of `addr`. Will be sorted to 
+* `S`: Tuple of the number of levels in each dimension, including the groundstate. Defaults
+    to a 1D spectrum with number of levels matching modes of `addr`. Will be sorted to
     make the first dimension the largest.
 * `η`: Define a custom aspect ratio for the trapping potential strengths, instead of deriving
-    from `S .- 1`. The values are always scaled relative to the first dimension, which sets 
+    from `S .- 1`. The values are always scaled relative to the first dimension, which sets
     the energy scale of the system, ``\\hbar\\omega_x``.
-* `g`: the (isotropic) interparticle interaction parameter. The value of `g` is assumed 
+* `g`: the (isotropic) interparticle interaction parameter. The value of `g` is assumed
     to be in trap units.
-* `interaction_only`: if set to `true` then the noninteracting single-particle terms are 
+* `interaction_only`: if set to `true` then the noninteracting single-particle terms are
     ignored. Useful if only energy shifts due to interactions are required.
 
 
@@ -113,9 +113,9 @@ struct HOCartesianEnergyConservedPerDim{
 end
 
 function HOCartesianEnergyConservedPerDim(
-        addr::BoseFS; 
+        addr::BoseFS;
         S = (num_modes(addr),),
-        η = nothing, 
+        η = nothing,
         g = 1.0,
         interaction_only = false
     )
@@ -147,7 +147,7 @@ function HOCartesianEnergyConservedPerDim(
     bigrange = 0:M-1
     # case n = M is the same as n = 0
     vmat = reshape(
-                [four_oscillator_integral_general(i-n, j+n, j, i; max_level = M-1) 
+                [four_oscillator_integral_general(i-n, j+n, j, i; max_level = M-1)
                     for i in bigrange, j in bigrange, n in bigrange],
                     M,M,M)
 
@@ -155,10 +155,11 @@ function HOCartesianEnergyConservedPerDim(
 end
 
 function Base.show(io::IO, h::HOCartesianEnergyConservedPerDim)
+    compact_addr = repr(h.addr, context=:compact => true) # compact print address
     flag = iszero(h.energies)
     # invert the scaling of u parameter
     g = h.u * 2 * sqrt(prod(h.aspect1))
-    print(io, "HOCartesianEnergyConservedPerDim($(h.addr); S=$(h.S), η=$(h.aspect1), g=$g, interaction_only=$flag)")
+    print(io, "HOCartesianEnergyConservedPerDim($(compact_addr); S=$(h.S), η=$(h.aspect1), g=$g, interaction_only=$flag)")
 end
 
 Base.:(==)(H::HOCartesianEnergyConservedPerDim, G::HOCartesianEnergyConservedPerDim) = all(map(p -> getproperty(H, p) == getproperty(G, p), propertynames(H)))
@@ -174,7 +175,7 @@ LOStructure(::Type{<:HOCartesianEnergyConservedPerDim}) = IsHermitian()
 function energy_transfer_diagonal(h::HOCartesianEnergyConservedPerDim{D}, omm::BoseOccupiedModeMap) where {D}
     result = 0.0
     states = CartesianIndices(h.S)    # 1-indexed
-    
+
     for i in eachindex(omm)
         mode_i, occ_i = omm[i].mode, omm[i].occnum
         idx_i = Tuple(states[mode_i])
@@ -188,7 +189,7 @@ function energy_transfer_diagonal(h::HOCartesianEnergyConservedPerDim{D}, omm::B
             idx_j = Tuple(states[mode_j])
             val = 4 * occ_i * occ_j
             result += prod(h.vtable[idx_i[d],idx_j[d],1] for d in 1:D) * val
-        end        
+        end
     end
     return result * h.u
 end
@@ -219,12 +220,12 @@ function num_offdiagonals(h::HOCartesianEnergyConservedPerDim, addr::BoseFS)
         p_i = omm[i]
         if p_i.occnum > 1
             _, valid_box_size  = largest_two_point_box(p_i.mode, p_i.mode, S)
-            noffs += valid_box_size 
+            noffs += valid_box_size
         end
         for j in 1:i-1
             p_j = omm[j]
             _, valid_box_size  = largest_two_point_box(p_i.mode, p_j.mode, S)
-            noffs += valid_box_size 
+            noffs += valid_box_size
         end
     end
     return noffs
@@ -240,9 +241,9 @@ Return the new address `new_add`, the prefactor `val`, the initial particle mode
 `mode_k` is implicit by energy conservation.
 """
 function energy_transfer_offdiagonal(
-        S::Tuple, 
-        addr::BoseFS, 
-        chosen::Int, 
+        S::Tuple,
+        addr::BoseFS,
+        chosen::Int,
         omm::BoseOccupiedModeMap = OccupiedModeMap(addr)
     )
     # find size of valid moves for each pair
@@ -267,14 +268,14 @@ function energy_transfer_offdiagonal(
 end
 
 function get_offdiagonal(
-        h::HOCartesianEnergyConservedPerDim{D,A}, 
-        addr::A, 
-        chosen::Int, 
+        h::HOCartesianEnergyConservedPerDim{D,A},
+        addr::A,
+        chosen::Int,
         omm::BoseOccupiedModeMap = OccupiedModeMap(addr)
     ) where {D,A}
 
     S = h.S
-    
+
     new_add, val, i, j, l = energy_transfer_offdiagonal(S, addr, chosen, omm)
 
     if val ≠ 0.0    # is this a safe check? maybe check if Δns == (0,...)
diff --git a/src/Hamiltonians/HubbardMom1DEP.jl b/src/Hamiltonians/HubbardMom1DEP.jl
index 61aadc848..2ba11098a 100644
--- a/src/Hamiltonians/HubbardMom1DEP.jl
+++ b/src/Hamiltonians/HubbardMom1DEP.jl
@@ -66,7 +66,7 @@ See also [`HubbardMom1D`](@ref), [`HubbardReal1DEP`](@ref),
 [`Transcorrelated1D`](@ref), [`Hamiltonians`](@ref).
 """
 struct HubbardMom1DEP{TT,M,AD<:AbstractFockAddress,U,T,V,D} <: AbstractHamiltonian{TT}
-    addr::AD # default starting address, should have M modes, U and T are model parameters
+    address::AD # default starting address, should have M modes, U and T are model parameters
     ks::SVector{M,TT} # values for k
     kes::SVector{M,TT} # values for kinetic energy
     ep::SVector{M,TT} # external potential
@@ -74,10 +74,10 @@ struct HubbardMom1DEP{TT,M,AD<:AbstractFockAddress,U,T,V,D} <: AbstractHamiltoni
 end
 
 function HubbardMom1DEP(
-    addr::Union{SingleComponentFockAddress,FermiFS2C};
+    address::Union{SingleComponentFockAddress,FermiFS2C};
     u=1.0, t=1.0, dispersion = hubbard_dispersion, v_ho=1.0,
 )
-    M = num_modes(addr)
+    M = num_modes(address)
     U, T, V = promote(float(u), float(t), float(v_ho))
     step = 2π/M
     if isodd(M)
@@ -91,22 +91,22 @@ function HubbardMom1DEP(
 
     potential = momentum_space_harmonic_potential(M, V)
 
-    return HubbardMom1DEP{typeof(U),M,typeof(addr),U,T,V,typeof(dispersion)}(
-        addr, ks, kes, potential, dispersion
+    return HubbardMom1DEP{typeof(U),M,typeof(address),U,T,V,typeof(dispersion)}(
+        address, ks, kes, potential, dispersion
     )
 end
 
 function starting_address(h::HubbardMom1DEP)
-    return h.addr
+    return h.address
 end
 
 dimension(::HubbardMom1DEP, address) = number_conserving_dimension(address)
 
-function get_offdiagonal(h::HubbardMom1DEP{<:Any,<:Any,F}, addr::F, i) where {F}
-    return offdiagonals(h, addr)[i]
+function get_offdiagonal(h::HubbardMom1DEP{<:Any,<:Any,F}, address::F, i) where {F}
+    return offdiagonals(h, address)[i]
 end
-function num_offdiagonals(h::HubbardMom1DEP{<:Any,<:Any,F}, addr::F) where {F}
-    return length(offdiagonals(h, addr))
+function num_offdiagonals(h::HubbardMom1DEP{<:Any,<:Any,F}, address::F) where {F}
+    return length(offdiagonals(h, address))
 end
 
 LOStructure(::Type{<:HubbardMom1DEP{<:Real}}) = IsHermitian()
@@ -115,20 +115,20 @@ Base.getproperty(h::HubbardMom1DEP, s::Symbol) = getproperty(h, Val(s))
 Base.getproperty(h::HubbardMom1DEP, ::Val{:ep}) = getfield(h, :ep)
 Base.getproperty(h::HubbardMom1DEP, ::Val{:ks}) = getfield(h, :ks)
 Base.getproperty(h::HubbardMom1DEP, ::Val{:kes}) = getfield(h, :kes)
-Base.getproperty(h::HubbardMom1DEP, ::Val{:addr}) = getfield(h, :addr)
+Base.getproperty(h::HubbardMom1DEP, ::Val{:address}) = getfield(h, :address)
 Base.getproperty(h::HubbardMom1DEP, ::Val{:dispersion}) = getfield(h, :dispersion)
 Base.getproperty(h::HubbardMom1DEP{<:Any,<:Any,<:Any,U}, ::Val{:u}) where {U} = U
 Base.getproperty(h::HubbardMom1DEP{<:Any,<:Any,<:Any,<:Any,T}, ::Val{:t}) where {T} = T
 Base.getproperty(h::HubbardMom1DEP{<:Any,<:Any,<:Any,<:Any,<:Any,V}, ::Val{:v_ho}) where {V} = V
 
 function Base.show(io::IO, h::HubbardMom1DEP)
-    compact_addr = repr(h.addr, context=:compact => true) # compact print address
+    compact_addr = repr(h.address, context=:compact => true) # compact print address
     print(io, "HubbardMom1DEP($compact_addr; ")
     print(io, "u=$(h.u), t=$(h.t), v_ho=$(h.v_ho), dispersion=$(h.dispersion))")
 end
 
 function Base.:(==)(a::HubbardMom1DEP, b::HubbardMom1DEP)
-    result = a.addr == b.addr && a.u == b.u && a.t == b.t && a.dispersion == b.dispersion
+    result = a.address == b.address && a.u == b.u && a.t == b.t && a.dispersion == b.dispersion
     return result && a.ep == b.ep && a.ks == b.ks && a.kes == b.kes
 end
 
@@ -145,13 +145,13 @@ end
     return h.u / 2M * momentum_transfer_diagonal(map_a, map_b)
 end
 
-@inline function diagonal_element(h::HubbardMom1DEP, addr::SingleComponentFockAddress)
-    map = OccupiedModeMap(addr)
+@inline function diagonal_element(h::HubbardMom1DEP, address::SingleComponentFockAddress)
+    map = OccupiedModeMap(address)
     return dot(h.kes, map) + momentum_transfer_diagonal(h, map) +
-        momentum_external_potential_diagonal(h.ep, addr, map)
+        momentum_external_potential_diagonal(h.ep, address, map)
 end
-@inline function diagonal_element(h::HubbardMom1DEP, addr::FermiFS2C)
-    c1, c2 = addr.components
+@inline function diagonal_element(h::HubbardMom1DEP, address::FermiFS2C)
+    c1, c2 = address.components
     map_a = OccupiedModeMap(c1)
     map_b = OccupiedModeMap(c2)
     return dot(h.kes, map_a) + dot(h.kes, map_b) +
diff --git a/src/Hamiltonians/HubbardReal1D.jl b/src/Hamiltonians/HubbardReal1D.jl
index 47c4cd243..8ce0f7266 100644
--- a/src/Hamiltonians/HubbardReal1D.jl
+++ b/src/Hamiltonians/HubbardReal1D.jl
@@ -29,6 +29,7 @@ function HubbardReal1D(addr; u=1.0, t=1.0)
 end
 
 function Base.show(io::IO, h::HubbardReal1D)
+    io = IOContext(io, :compact => true)
     print(io, "HubbardReal1D(")
     show(io, h.add)
     print(io, "; u=$(h.u), t=$(h.t))")
diff --git a/src/Hamiltonians/HubbardReal1DEP.jl b/src/Hamiltonians/HubbardReal1DEP.jl
index 45fa6a38e..8b7ae7c67 100644
--- a/src/Hamiltonians/HubbardReal1DEP.jl
+++ b/src/Hamiltonians/HubbardReal1DEP.jl
@@ -45,21 +45,22 @@ lattice.
 
 """
 struct HubbardReal1DEP{TT,A<:AbstractFockAddress,U,T,M} <: AbstractHamiltonian{TT}
-    add::A
+    address::A
     ep::SVector{M,TT}
 end
 
-function HubbardReal1DEP(addr::SingleComponentFockAddress{<:Any,M}; u=1.0, t=1.0, v_ho=1.0) where M
+function HubbardReal1DEP(address::SingleComponentFockAddress{<:Any,M}; u=1.0, t=1.0, v_ho=1.0) where M
     U, T, V = promote(float(u), float(t), float(v_ho))
     # js = range(1-cld(M,2); length=M)
     is = range(-fld(M,2); length=M) # [-M÷2, M÷2) including left boundary
     js = shift_lattice(is) # shifted such that js[1] = 0
     potential = SVector{M}(V*j^2 for j in js)
-    return HubbardReal1DEP{typeof(U),typeof(addr),U,T,M}(addr, potential)
+    return HubbardReal1DEP{typeof(U),typeof(address),U,T,M}(address, potential)
 end
 
 function Base.show(io::IO, h::HubbardReal1DEP)
-    print(io, "HubbardReal1DEP($(h.add); u=$(h.u), t=$(h.t), v_ho=$(h.ep[2]))")
+    compact_addr = repr(h.address, context=:compact => true) # compact print address
+    print(io, "HubbardReal1DEP($(compact_addr); u=$(h.u), t=$(h.t), v_ho=$(h.ep[2]))")
 end
 
 LOStructure(::Type{<:HubbardReal1DEP{<:Real}}) = IsHermitian()
@@ -67,10 +68,10 @@ LOStructure(::Type{<:HubbardReal1DEP{<:Real}}) = IsHermitian()
 Base.getproperty(h::HubbardReal1DEP, s::Symbol) = getproperty(h, Val(s))
 Base.getproperty(h::HubbardReal1DEP{<:Any,<:Any,U}, ::Val{:u}) where U = U
 Base.getproperty(h::HubbardReal1DEP{<:Any,<:Any,<:Any,T}, ::Val{:t}) where T = T
-Base.getproperty(h::HubbardReal1DEP, ::Val{:add}) = getfield(h, :add)
+Base.getproperty(h::HubbardReal1DEP, ::Val{:address}) = getfield(h, :address)
 Base.getproperty(h::HubbardReal1DEP, ::Val{:ep}) = getfield(h, :ep)
 
-starting_address(h::HubbardReal1DEP) = h.add
+starting_address(h::HubbardReal1DEP) = h.address
 
 dimension(::HubbardReal1DEP, address) = number_conserving_dimension(address)
 
@@ -85,7 +86,7 @@ function diagonal_element(h::HubbardReal1DEP, address::SingleComponentFockAddres
     end
 end
 
-function get_offdiagonal(h::HubbardReal1DEP, add::SingleComponentFockAddress, chosen)
-    naddress, onproduct = hopnextneighbour(add, chosen)
+function get_offdiagonal(h::HubbardReal1DEP, address::SingleComponentFockAddress, chosen)
+    naddress, onproduct = hopnextneighbour(address, chosen)
     return naddress, - h.t * onproduct
 end
diff --git a/src/Hamiltonians/HubbardRealSpace.jl b/src/Hamiltonians/HubbardRealSpace.jl
index 9e25c64a6..dd1a75fb8 100644
--- a/src/Hamiltonians/HubbardRealSpace.jl
+++ b/src/Hamiltonians/HubbardRealSpace.jl
@@ -255,6 +255,7 @@ warn_fermi_interaction(_, _) = nothing
 LOStructure(::Type{<:HubbardRealSpace}) = IsHermitian()
 
 function Base.show(io::IO, h::HubbardRealSpace{C}) where C
+    io = IOContext(io, :compact => true)
     println(io, "HubbardRealSpace(")
     println(io, "  ", starting_address(h), ",")
     println(io, "  geometry = ", h.geometry, ",")
diff --git a/src/Hamiltonians/ParitySymmetry.jl b/src/Hamiltonians/ParitySymmetry.jl
index f5f72b292..d4267a50f 100644
--- a/src/Hamiltonians/ParitySymmetry.jl
+++ b/src/Hamiltonians/ParitySymmetry.jl
@@ -24,7 +24,7 @@ parity, respectively, where `ᾱ == reverse(α)`.
 
 ```jldoctest
 julia> ham = HubbardReal1D(BoseFS(0,2,1))
-HubbardReal1D(BoseFS{3,3}(0, 2, 1); u=1.0, t=1.0)
+HubbardReal1D(fs"|0 2 1⟩"; u=1.0, t=1.0)
 
 julia> size(Matrix(ham))
 (10, 10)
@@ -116,4 +116,3 @@ end
 function diagonal_element(h::ParitySymmetry, add)
     return diagonal_element(h.hamiltonian, add)
 end
-
diff --git a/src/Hamiltonians/Transcorrelated1D.jl b/src/Hamiltonians/Transcorrelated1D.jl
index 0187f1f54..514b3eead 100644
--- a/src/Hamiltonians/Transcorrelated1D.jl
+++ b/src/Hamiltonians/Transcorrelated1D.jl
@@ -89,6 +89,7 @@ function Transcorrelated1D(address; t=1.0, v=1.0, v_ho=0.0, cutoff=1, three_body
 end
 
 function Base.show(io::IO, h::Transcorrelated1D)
+    io = IOContext(io, :compact => true)
     print(io, "Transcorrelated1D(", starting_address(h), ", t=$(h.t), v=$(h.v)")
     if !iszero(h.v_ho)
         print(io, ", v_ho=$(h.v_ho))")
diff --git a/src/lomc.jl b/src/lomc.jl
index 5bd963bdb..a1a8d4f22 100644
--- a/src/lomc.jl
+++ b/src/lomc.jl
@@ -51,9 +51,9 @@ Some metadata is automatically added to the report `df` including
 # Example
 
 ```jldoctest
-julia> add = BoseFS(1,2,3);
+julia> address = BoseFS(1,2,3);
 
-julia> hamiltonian = HubbardReal1D(add);
+julia> hamiltonian = HubbardReal1D(address);
 
 julia> df1, state = lomc!(hamiltonian; targetwalkers=500, laststep=100);
 
@@ -69,7 +69,7 @@ julia> using DataFrames; metadata(df2, "info") # retrieve custom metadata
 "cont"
 
 julia> metadata(df2, "hamiltonian") # some metadata is automatically added
-"HubbardReal1D(BoseFS{6,3}(1, 2, 3); u=1.0, t=1.0)"
+"HubbardReal1D(fs\\"|1 2 3⟩\\"; u=1.0, t=1.0)"
 ```
 
 # Further keyword arguments and defaults: