Module `Hamiltonians`

This module contains definitions of Hamiltonians, in particular specific physical models of interest. These are organised by means of an interface around the abstract type AbstractHamiltonian, in the spirit of the AbstractArray interface as discussed in the Julia Documentation.

Rimu.Hamiltonians — Module

module Hamiltonians

This module defines Hamiltonian types and functions for working with Hamiltonians.

Exported concrete Hamiltonian types

Real space Hubbard models

Momentum space Hubbard models

Harmonic oscillator models

Other

Wrappers

Observables

Interface for working with Hamiltonians

AbstractHamiltonian: defined in the module Interfaces

source

Usage with FCIQMC and exact diagonalisation

In order to define a specific model Hamiltonian with relevant parameters for the model, instantiate the model like this in the input file:

hubb = HubbardReal1D(BoseFS((1,2,0,3)); u=1.0, t=1.0)

The Hamiltonian hubb is now ready to be used for FCIQMC in lomc! and for exact diagonalisation with KrylovKit.jl directly, or after transforming into a sparse matrix first with

using SparseArrays
sh = sparse(hubb)

or into a full matrix with

using LinearAlgebra
fh = Matrix(hubb)

This functionality relies on

Rimu.Hamiltonians.BasisSetRep — Type

BasisSetRep(
    h::AbstractHamiltonian, addr=starting_address(h);
    sizelim=10^6, nnzs, cutoff, filter, sort=false, kwargs...
)
BasisSetRep(h::AbstractHamiltonian, addresses::AbstractVector; kwargs...)

Eagerly construct the basis set representation of the operator h with all addresses reachable from addr. Instead of a single address, a vector of addresses can be passed.

An ArgumentError is thrown if dimension(h) > sizelim in order to prevent memory overflow. Set sizelim = Inf in order to disable this behaviour.

Providing the number nnzs of expected calculated matrix elements and col_hint for the estimated number of nonzero off-diagonal matrix elements in each matrix column may improve performance.

Providing an energy cutoff will skip the columns and rows with diagonal elements greater than cutoff. Alternatively, an arbitrary filter function can be used instead. Addresses passed as arguments are not filtered. To generate the matrix truncated to the subspace spanned by the addresses, use filter = Returns(false).

Setting sort to true will sort the matrix rows and columns. This is useful when the order of the columns matters, e.g. when comparing matrices. Any additional keyword arguments are passed on to Base.sortperm.

Fields

sm: sparse matrix representing h in the basis basis
basis: vector of addresses
h: the Hamiltonian

Example

julia> h = HubbardReal1D(BoseFS(1,0,0));

julia> bsr = BasisSetRep(h)
BasisSetRep(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:
  0.0  -1.0  -1.0
 -1.0   0.0  -1.0
 -1.0  -1.0   0.0

julia> BasisSetRep(h, bsr.basis[1:2]; filter = Returns(false)) # passing addresses and truncating
BasisSetRep(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:
  0.0  -1.0
 -1.0   0.0

julia> using LinearAlgebra; eigvals(Matrix(bsr)) # eigenvalues
3-element Vector{Float64}:
 -1.9999999999999996
  0.9999999999999997
  1.0000000000000002

julia> ev = eigvecs(Matrix(bsr))[:,1] # ground state eigenvector
3-element Vector{Float64}:
 -0.5773502691896257
 -0.5773502691896255
 -0.5773502691896257

julia> DVec(zip(bsr.basis,ev)) # ground state as DVec
DVec{BoseFS{1, 3, BitString{3, 1, UInt8}},Float64} with 3 entries, style = IsDeterministic{Float64}()
  fs"|0 0 1⟩" => -0.57735
  fs"|0 1 0⟩" => -0.57735
  fs"|1 0 0⟩" => -0.57735

Has methods for dimension, sparse, Matrix, starting_address.

source

SparseArrays.sparse — Function

sparse(h::AbstractHamiltonian, addr=starting_address(h); kwargs...)
sparse(bsr::BasisSetRep)

Return a sparse matrix representation of h or bsr. kwargs are passed to BasisSetRep.

See BasisSetRep.

source

Base.Matrix — Type

Matrix(h::AbstractHamiltonian, addr=starting_address(h); sizelim=10^4, kwargs...)
Matrix(bsr::BasisSetRep)

Return a dense matrix representation of h or bsr. kwargs are passed to BasisSetRep.

See BasisSetRep.

source

If only the basis is required and not the matrix representation it is more efficient to use

Rimu.Hamiltonians.build_basis — Function

build_basis(
    ham, address=starting_address(ham);
    cutoff, filter, sizelim, sort=false, kwargs...
) -> basis
build_basis(ham, addresses::AbstractVector; kwargs...)

Get all basis element of a linear operator ham that are reachable (via non-zero matrix elements) from the address address, returned as a vector. Instead of a single address, a vector of addresses can be passed. Does not return the matrix, for that purpose use BasisSetRep.

Providing an energy cutoff will skip addresses with diagonal elements greater than cutoff. Alternatively, an arbitrary filter function can be used instead. Addresses passed as arguments are not filtered. A maximum basis size sizelim can be set which will throw an error if the expected dimension of ham is larger than sizelim. This may be useful when memory may be a concern. These options are disabled by default.

Setting sort to true will sort the basis. Any additional keyword arguments are passed on to Base.sort!.

source

Model Hamiltonians

Here is a list of fully implemented model Hamiltonians. There are several variants of the Hubbard model in real and momentum space, as well as some other models.

Real space Hubbard models

Rimu.Hamiltonians.HubbardReal1D — Type

HubbardReal1D(address; u=1.0, t=1.0)

Implements a one-dimensional Bose Hubbard chain in real space.

\[\hat{H} = -t \sum_{\langle i,j\rangle} a_i^† a_j + \frac{u}{2}\sum_i n_i (n_i-1)\]

Arguments

address: the starting address, defines number of particles and sites.
u: the interaction parameter.
t: the hopping strength.

See also

source

Rimu.Hamiltonians.BoseHubbardReal1D2C — Type

BoseHubbardReal1D2C(address::BoseFS2C; ua=1.0, ub=1.0, ta=1.0, tb=1.0, v=1.0)

Implements a two-component one-dimensional Bose Hubbard chain in real space.

\[\hat{H} = \hat{H}_a + \hat{H}_b + V\sum_{i} n_{a_i}n_{b_i}\]

Arguments

address: the starting address, defines number of particles and sites.
ua: the on-site interaction parameter parameter for Hamiltonian a.
ub: the on-site interaction parameter parameter for Hamiltonian b.
ta: the hopping strength for Hamiltonian a.
tb: the hopping strength for Hamiltonian b.
v: the inter-species interaction parameter V.

See also

source

Rimu.Hamiltonians.HubbardReal1DEP — Type

HubbardReal1DEP(address; u=1.0, t=1.0, v_ho=1.0)

Implements a one-dimensional Bose Hubbard chain in real space with external potential.

\[\hat{H} = -t \sum_{\langle i,j\rangle} a_i^† a_j + \sum_i ϵ_i n_i + \frac{u}{2}\sum_i n_i (n_i-1)\]

Arguments

address: the starting address, defines number of particles and sites.
u: the interaction parameter.
t: the hopping strength.
v_ho: strength of the external harmonic oscillator potential $ϵ_i = v_{ho} i^2$.

The first index is i=0 and the maximum of the potential occurs in the centre of the lattice.

See also

source

Rimu.Hamiltonians.HubbardRealSpace — Type

HubbardRealSpace(address; geometry=PeriodicBoundaries(M,), t=ones(C), u=ones(C, C), v=zeros(C, D))

Hubbard model in real space. Supports single or multi-component Fock state addresses (with C components) and various (rectangular) lattice geometries in D dimensions.

\[ \hat{H} = -\sum_{\langle i,j\rangle,σ} t_σ a^†_{iσ} a_{jσ} + \frac{1}{2}\sum_{i,σ} u_{σσ} n_{iσ} (n_{iσ} - 1) + \sum_{i,σ≠τ}u_{στ} n_{iσ} n_{iτ}\]

If v is nonzero then this calculates $\hat{H} + \hat{V}$ by adding the harmonic trapping potential

\[ \hat{V} = \sum_{i,σ,d} v_{σd} x_{di}^2 n_{iσ}\]

where $x_{di}$ is the distance of site $i$ from the centre of the trap along dimension $d$.

Address types

BoseFS: Single-component Bose-Hubbard model.
FermiFS: Single-component Fermi-Hubbard model.
CompositeFS: For multi-component models.

Note that a single component of fermions cannot interact with itself. A warning is produced if addressis incompatible with the interaction parameters u.

Geometries

Implemented LatticeGeometrys for keyword geometry

Default is geometry=PeriodicBoundaries(M,), i.e. a one-dimensional lattice with the number of sites M inferred from the number of modes in address.

Other parameters

t: the hopping strengths. Must be a vector of length C. The i-th element of the vector corresponds to the hopping strength of the i-th component.
u: the on-site interaction parameters. Must be a symmetric matrix. u[i, j] corresponds to the interaction between the i-th and j-th component. u[i, i] corresponds to the interaction of a component with itself. Note that u[i,i] must be zero for fermionic components.
v: the trap potential strengths. Must be a matrix of size C × D. v[i,j] is the strength of the trap for component i in the jth dimension.

source

Rimu.Hamiltonians.ExtendedHubbardReal1D — Type

ExtendedHubbardReal1D(address; u=1.0, v=1.0, t=1.0)

Implements the extended Hubbard model on a one-dimensional chain in real space.

\[\hat{H} = -t \sum_{\langle i,j\rangle} a_i^† a_j + \frac{u}{2}\sum_i n_i (n_i-1) + v \sum_{\langle i,j\rangle} n_i n_j\]

Arguments

address: the starting address.
u: on-site interaction parameter
v: the next-neighbor interaction
t: the hopping strength

source

Momentum space Hubbard models

Rimu.Hamiltonians.HubbardMom1D — Type

HubbardMom1D(address; u=1.0, t=1.0, dispersion=hubbard_dispersion)

Implements a one-dimensional Bose Hubbard chain in momentum space.

\[\hat{H} = \sum_{k} ϵ_k n_k + \frac{u}{M}\sum_{kpqr} a^†_{r} a^†_{q} a_p a_k δ_{r+q,p+k}\]

Arguments

address: the starting address, defines number of particles and sites.
u: the interaction parameter.
t: the hopping strength.
dispersion: defines $ϵ_k =$t*dispersion(k)
- hubbard_dispersion: $ϵ_k = -2t \cos(k)$
- continuum_dispersion: $ϵ_k = tk^2$

See also

source

Rimu.Hamiltonians.BoseHubbardMom1D2C — Type

BoseHubbardMom1D2C(add::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.

\[\hat{H} = \hat{H}_a + \hat{H}_b + \frac{V}{M}\sum_{kpqr} b^†_{r} a^†_{q} b_p a_k δ_{r+q,p+k}\]

Arguments

add: the starting address.
ua: the u parameter for Hamiltonian a.
ub: the u parameter for Hamiltonian b.
ta: the t parameter for Hamiltonian a.
tb: the t parameter for Hamiltonian b.
v: the inter-species interaction parameter V.

Further keyword arguments are passed on to the constructor of HubbardMom1D.

See also

source

Rimu.Hamiltonians.HubbardMom1DEP — Type

HubbardMom1DEP(address; u=1.0, t=1.0, v_ho=1.0, dispersion=hubbard_dispersion)

Implements a one-dimensional Bose Hubbard chain in momentum space with harmonic external potential.

\[Ĥ = \sum_{k} ϵ_k n_k + \frac{u}{M}\sum_{kpqr} a^†_{r} a^†_{q} a_p a_k δ_{r+q,p+k} + V̂_\mathrm{ho} ,\]

where

\[\begin{aligned} V̂_\mathrm{ho} & = \frac{1}{M} \sum_{p,q} \mathrm{DFT}[V_{ext}]_{p-q} \, a^†_{p} a_q ,\\ V_\mathrm{ext}(x) &= v_\mathrm{ho} \,x^2 , \end{aligned}\]

is an external harmonic potential in momentum space, $\mathrm{DFT}[…]_k$ is a discrete Fourier transform performed by fft()[k%M + 1], and M == num_modes(address).

Arguments

address: the starting address, defines number of particles and sites.
u: the interaction parameter.
t: the hopping strength.
dispersion: defines $ϵ_k =$t*dispersion(k)
- hubbard_dispersion: $ϵ_k = -2t \cos(k)$
- continuum_dispersion: $ϵ_k = tk^2$
v_ho: strength of the external harmonic oscillator potential $v_\mathrm{ho}$.

source

Harmonic oscillator models

Rimu.Hamiltonians.HOCartesianContactInteractions — Type

HOCartesianContactInteractions(addr; S, η, g = 1.0, interaction_only = false, block_by_level = true)

Implements a bosonic harmonic oscillator in Cartesian basis with contact interactions

\[\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

\[ \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.

By default the offdiagonal elements due to the interactions are consistent with first-order degenerate perturbation theory

\[ V_{\mathbf{ijkl}} = \delta_{\epsilon_\mathbf{i} + \epsilon_\mathbf{j}} ^{\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 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 fewer but larger blocks (the blocks are still distinguished by parity of basis states). Alternatively, see HOCartesianEnergyConservedPerDim 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 basis functions, implemented in four_oscillator_integral_general.

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 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 the energy scale of the system, $\hbar\omega_x$.
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 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 comparing their noninteracting energy.

Warning

num_offdiagonals is a bad estimate for this Hamiltonian. Take care when building a matrix or using QMC methods. Use get_all_blocks first then pass option col_hint = block_size to BasisSetRep to safely build the matrix.

source

Rimu.Hamiltonians.HOCartesianEnergyConservedPerDim — Type

HOCartesianEnergyConservedPerDim(addr; S, η, g = 1.0, interaction_only = false)

Implements a bosonic harmonic oscillator in Cartesian basis with contact interactions

\[\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 for a model that conserves total energy.

\[ \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{ijkl}}$ are for a contact interaction calculated in this basis using first-order degenerate perturbation theory.

\[ 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, 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 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 the energy scale of the system, $\hbar\omega_x$.
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 ignored. Useful if only energy shifts due to interactions are required.

source

Rimu.Hamiltonians.HOCartesianCentralImpurity — Type

HOCartesianCentralImpurity(addr; kwargs...)

Hamiltonian of non-interacting particles in an arbitrary harmonic trap with a delta-function potential at the centre, with strength g,

\[\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}.\]

\[ \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

\[ 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 Ints.
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.

source

Other

Rimu.Hamiltonians.MatrixHamiltonian — Type

MatrixHamiltonian(
    mat::AbstractMatrix{T};
    starting_address::Int = starting_address(mat)
) <: AbstractHamiltonian{T}

Wrap an abstract matrix mat as an AbstractHamiltonian object. Works with stochastic methods of lomc!() and DVec. Optionally, a valid index can be provided as the starting_address.

Specialised methods are implemented for sparse matrices of type AbstractSparseMatrixCSC. One based indexing is required for the matrix mat.

source

Rimu.Hamiltonians.Transcorrelated1D — Type

Transcorrelated1D(address; t=1.0, v=1.0, v_ho=0.0, cutoff=1, three_body_term=true)

Implements a transcorrelated Hamiltonian for contact interactions in one dimensional momentum space from Jeszenski et al. (2018). Currently limited to two component fermionic addresses.

\[\begin{aligned} \tilde{H} &= t \sum_{kσ}k^2 n_{k,σ} \\ &\quad + \sum_{pqkσσ'} T_{pqk} a^†_{p-k,σ} a^†_{q+k,σ'} a_{q,σ'} a_{p,σ} \\ &\quad + \sum_{pqskk'σσ'} Q_{kk'}a^†_{p-k,σ} a^†_{q+k,σ} a^†_{s+k-k',σ'} a_{s,σ'} a_{q,σ} a_{p,σ} \\ &\quad + V̂_\mathrm{ho} \end{aligned}\]

where

\[\begin{aligned} \tilde{u}(k) &= \begin{cases} -\frac{2}{k^2} &\mathrm{if\ } |k| ≥ k_c\\ 0 & \mathrm{otherwise} \end{cases} \\ T_{pqk} &= \frac{v}{M} + \frac{2v}{M}\left[k^2\tilde{u}(k) - (p - q)k\tilde{u}(k)\right] + \frac{2v^2}{t}W(k)\\ W(k) &= \frac{1}{M^2}\sum_{q} (k - q)q\, \tilde{u}(q)\,\tilde{u}(k - q) \\ Q_{kl} &= -\frac{v^2}{t M^2}k \tilde{u}(k)\,l\tilde{u}(l), \end{aligned}\]

Arguments

address: The starting address, defines number of particles and sites.
v: The interaction parameter.
t: The kinetic energy prefactor.
v_ho: Strength of the external harmonic oscillator potential $V̂_\mathrm{ho}$. See HubbardMom1DEP.
cutoff controls $k_c$ in equations above. Note: skipping generating off-diagonal elements below the cutoff is not implemented - zero-valued elements are returned instead.
three_body_term: If set to false, generating three body excitations is skipped. Note: when disabling three body terms, cutoff should be set to a higher value for best results.

See also

source

Convenience functions

Rimu.Hamiltonians.rayleigh_quotient — Function

rayleigh_quotient(H, v)

\[\frac{⟨ v | H | v ⟩}{⟨ v|v ⟩}\]

source

Rimu.Hamiltonians.momentum — Function

momentum(ham::AbstractHamiltonian)

Momentum as a linear operator in Fock space. Pass a Hamiltonian ham in order to convey information about the Fock basis.

Note: momentum is currently only defined on HubbardMom1D.

Example

julia> add = BoseFS((1, 0, 2, 1, 2, 1, 1, 3));


julia> ham = HubbardMom1D(add; u = 2.0, t = 1.0);


julia> mom = momentum(ham);


julia> diagonal_element(mom, add) # calculate the momentum of a single configuration
-1.5707963267948966

julia> v = DVec(add => 10; capacity=1000);


julia> rayleigh_quotient(mom, v) # momentum expectation value for state vector `v`
-1.5707963267948966

source

Rimu.Hamiltonians.hubbard_dispersion — Function

hubbard_dispersion(k)

Dispersion relation for HubbardMom1D. Returns -2cos(k).

Hamiltonian wrappers

The following Hamiltonians are constructed from an existing Hamiltonian instance and change its behaviour:

Rimu.Hamiltonians.GutzwillerSampling — Type

GutzwillerSampling(::AbstractHamiltonian; g)

Wrapper over any AbstractHamiltonian that implements Gutzwiller sampling. In this importance sampling scheme the Hamiltonian is modified as follows

\[\tilde{H}_{ij} = H_{ij} e^{-g(H_{ii} - H_{jj})} .\]

This way off-diagonal spawns to higher-energy configurations are discouraged and spawns to lower-energy configurations encouraged for positive g.

Constructor

GutzwillerSampling(::AbstractHamiltonian, g)
GutzwillerSampling(::AbstractHamiltonian; g)

After construction, we can access the underlying Hamiltonian with G.hamiltonian and the g parameter with G.g.

Example

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)

julia> G = GutzwillerSampling(H, g=0.3)
GutzwillerSampling(HubbardMom1D(BoseFS{3,3}(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)

julia> get_offdiagonal(G, BoseFS(2, 1, 0), 1)
(BoseFS{3,3}(1, 0, 2), 0.8131393194811987)

Observables

To calculate observables, pass the transformed Hamiltonian G to AllOverlaps with keyword argument transform=G.

source

Rimu.Hamiltonians.GuidingVectorSampling — Type

GuidingVectorSampling

Wrapper over any AbstractHamiltonian that implements guided vector a.k.a. guided wave function sampling. In this importance sampling scheme the Hamiltonian is modified as follows.

\[\tilde{H}_{ij} = v_i H_{ij} v_j^{-1}\]

and where v is the guiding vector. v_i and v_j are also thresholded to avoid dividing by zero (see below).

Constructors

GuidingVectorSampling(::AbstractHamiltonian, vector, eps)
GuidingVectorSampling(::AbstractHamiltonian; vector, eps)

eps is a thresholding parameter used to avoid dividing by zero; all values below eps are set to eps. It is recommended that eps is in the same value range as the guiding vector. The default value is set to eps=norm(v, Inf) * 1e-2

After construction, we can access the underlying hamiltonian with G.hamiltonian, the eps parameter with G.eps, and the guiding vector with G.vector.

Example

julia> H = HubbardReal1D(BoseFS(1,1,1); u=6.0, t=1.0);

julia> v = DVec(starting_address(H) => 10; capacity=1);

julia> G = GuidingVectorSampling(H, v, 0.1);

julia> get_offdiagonal(H, starting_address(H), 4)
(BoseFS{3,3}(2, 0, 1), -1.4142135623730951)

julia> get_offdiagonal(G, starting_address(G), 4)
(BoseFS{3,3}(2, 0, 1), -0.014142135623730952)

Observables

To calculate observables, pass the transformed Hamiltonian G to AllOverlaps with keyword argument transform=G.

source

Rimu.Hamiltonians.ParitySymmetry — Type

ParitySymmetry(ham::AbstractHamiltonian{T}; even=true) <: AbstractHamiltonian{T}

Impose even or odd parity on all states and the Hamiltonian ham as controlled by the keyword argument even. Parity symmetry of the Hamiltonian is assumed. For some Hamiltonians, ParitySymmetry reduces the size of the Hilbert space by half.

ParitySymmetry performs a unitary transformation, leaving the eigenvalues unchanged and preserving the LOStructure. This is achieved by changing the basis set to states with defined parity. Effectively, a non-even address $|α⟩$ is replaced by $\frac{1}{√2}(|α⟩ ± |ᾱ⟩)$ for even and odd parity, respectively, where ᾱ == reverse(α).

Notes

This modifier currently only works on starting_addresss with an odd number of modes.
For odd parity, the starting_address of the underlying Hamiltonian cannot be symmetric.
If parity is not a symmetry of the Hamiltonian ham then the result is undefined.
ParitySymmetry works by modifying the offdiagonals iterator.

julia> ham = HubbardReal1D(BoseFS(0,2,1))
HubbardReal1D(BoseFS{3,3}(0, 2, 1); u=1.0, t=1.0)

julia> size(Matrix(ham))
(10, 10)

julia> size(Matrix(ParitySymmetry(ham)))
(6, 6)

julia> size(Matrix(ParitySymmetry(ham; odd=true)))
(4, 4)

julia> eigvals(Matrix(ham))[1] ≈ eigvals(Matrix(ParitySymmetry(ham)))[1]
true

Observables

Observables are AbstractHamiltonians that represent a physical observable. Their ground state expectation values can be sampled by passing them into AllOverlaps.

Rimu.Hamiltonians.G2MomCorrelator — Type

G2MomCorrelator(d::Int,c=:cross) <: AbstractHamiltonian{ComplexF64}

Two-body correlation operator representing the density-density correlation at distance d of a two component system in a momentum-space Fock-state basis. It returns a Complex value.

Correlation across two components:

\[\hat{G}^{(2)}(d) = \frac{1}{M}\sum_{spqr=1}^M e^{-id(p-q)2π/M} a^†_{s} b^†_{p} b_q a_r δ_{s+p,q+r}\]

Correlation within a single component:

\[\hat{G}^{(2)}(d) = \frac{1}{M}\sum_{spqr=1}^M e^{-id(p-q)2π/M} a^†_{s} a^†_{p} a_q a_r δ_{s+p,q+r}\]

The diagonal element, where (p-q)=0, is

\[\frac{1}{M}\sum_{k,p=1}^M a^†_{k} b^†_{p} b_p a_k .\]

Arguments

d::Integer: the distance between two particles.
c: possible instructions: :cross: default instruction, computing correlation between particles across two components; :first: computing correlation between particles within the first component; :second: computing correlation between particles within the second component. These are the only defined instructions, using anything else will produce errors.

To use on a one-component system

For a system with only one component, e.g. with BoseFS, the second argument c is irrelevant and can be any of the above instructions, one could simply skip this argument and let it be the default value.

See also

source

Rimu.Hamiltonians.G2RealCorrelator — Type

G2RealCorrelator(d::Int) <: AbstractHamiltonian{Float64}

Two-body operator for density-density correlation between sites separated by d with 0 ≤ d < M.

\[ \hat{G}^{(2)}(d) = \frac{1}{M} \sum_i^M \hat{n}_i (\hat{n}_{i+d} - \delta_{0d}).\]

Assumes a one-dimensional lattice with periodic boundary conditions where

\[ \hat{G}^{(2)}(-M/2 \leq d < 0) = \hat{G}^{(2)}(|d|),\]

\[ \hat{G}^{(2)}(M/2 < d < M) = \hat{G}^{(2)}(M - d),\]

and normalisation

\[ \sum_{d=0}^{M-1} \langle \hat{G}^{(2)}(d) \rangle = \frac{N (N-1)}{M}.\]

For multicomponent basis, calculates correlations between all particles equally, equivalent to stacking all components into a single Fock state.

Arguments

d::Integer: distance between sites.

See also

source

Rimu.Hamiltonians.SuperfluidCorrelator — Type

SuperfluidCorrelator(d::Int) <: AbstractHamiltonian{Float64}

Operator for extracting superfluid correlation between sites separated by a distance d with 0 ≤ d < M:

\[ \hat{C}_{\text{SF}}(d) = \frac{1}{M} \sum_{i}^{M} a_{i}^{\dagger} a_{i + d}\]

Assumes a one-dimensional lattice with $M$ sites and periodic boundary conditions. $M$ is also the number of modes in the Fock state address.

Usage

Superfluid correlations can be extracted from a Monte Carlo calculation by wrapping SuperfluidCorrelator with AllOverlaps and passing into lomc! with the replica keyword argument. For an example with a similar use of G2RealCorrelator see G2 Correlator Example.

source

Rimu.Hamiltonians.StringCorrelator — Type

StringCorrelator(d::Int) <: AbstractHamiltonian{Float64}

Operator for extracting string correlation between lattice sites on a one-dimensional Hubbard lattice separated by a distance d with 0 ≤ d < M

\[ \hat{C}_{\text{string}}(d) = \frac{1}{M} \sum_{j}^{M} \delta n_j (e^{i \pi \sum_{j \leq k < j + d} \delta n_k}) \delta n_{j+d}\]

Here, $\delta \hat{n}_j = \hat{n}_j - \bar{n}$ is the boson number deviation from the mean filling number and $\bar{n} = N/M$ is the mean filling number of lattice sites with $N$ particles and $M$ lattice sites (or modes).

Assumes a one-dimensional lattice with periodic boundary conditions. For usage see SuperfluidCorrelator and AllOverlaps.

source

Rimu.Hamiltonians.DensityMatrixDiagonal — Type

DensityMatrixDiagonal(mode; component=0) <: AbstractHamiltonian

Represent a diagonal element of the single-particle density:

\[\hat{n}_{i,σ} = \hat a^†_{i,σ} \hat a_{i,σ}\]

where $i$ is the mode and $σ$ is the component. If component is zero, the sum over all components is computed.

See also

source

Rimu.Hamiltonians.Momentum — Type

Momentum(component=0; fold=true) <: AbstractHamiltonian

The momentum operator $\hat{p}$.

The component argument controls which component of the address is taken into consideration. A value of 0 sums the contributions of all components. If fold is true, the momentum is folded into the Brillouin zone.

julia> add = BoseFS((1, 0, 2, 1, 2, 1, 1, 3))
BoseFS{11,8}(1, 0, 2, 1, 2, 1, 1, 3)

julia> v = DVec(add => 10);

julia> rayleigh_quotient(Momentum(), DVec(add => 1))
-2.0

julia> rayleigh_quotient(Momentum(fold=false), DVec(add => 1))
14.0

source

Rimu.Hamiltonians.AxialAngularMomentumHO — Type

AxialAngularMomentumHO(S; z_dim = 3, addr = BoseFS(prod(S))) <: AbstractHamiltonian

Angular momentum operator for application to Cartesian harmonic oscillator basis, see HOCartesianContactInteractions or HOCartesianEnergyConservedPerDim. Represents the projection of angular momentum onto z-axis:

\[\hat{L}_z = i \hbar \sum_{j=1}^N \left( b_x b_y^\dag - b_y b_x^\dag \right),\]

where $b_x^\dag$ and $b_x$ are raising and lowering (ladder) operators for a harmonic oscillator in the $x$ dimension, and simlarly for $y$.

This is implemented for an $N$ particle Fock space with creation and annihilation operators as

\[\frac{1}{\hbar} \hat{L}_z = i \sum_{n_x=1}^{M_x} \sum_{n_y=1}^{M_y} \left( a_{n_x-1,n_y+1}^\dag - a_{n_x+1,n_y-1}^\dag \right) a_{n_x, n_y}.\]

in units of $\hbar$.

Argument S is a tuple defining the range of Cartesian modes in each dimension and their mapping to Fock space modes in a SingleComponentFockAddress. If S indicates a 3D system the z dimension can be changed by setting z_dim; S should be be isotropic in the remaining x-y plane, i.e. must have S[x_dim] == S[y_dim]. The starting address addr only needs to satisfy num_modes(addr) == prod(S).

Example

Calculate the overlap of two Fock addresses interpreted as harmonic oscillator states in a 2D Cartesian basis

julia> S = (2,2)
(2, 2)

julia> Lz = AxialAngularMomentumHO(S)
AxialAngularMomentumHO((2, 2); z_dim = 3, addr = BoseFS{0,4}(0, 0, 0, 0))

julia> v = DVec(BoseFS(prod(S), 2 => 1) => 1.0)
DVec{BoseFS{1, 4, BitString{4, 1, UInt8}},Float64} with 1 entry, style = IsDeterministic{Float64}()
  fs"|0 1 0 0⟩" => 1.0

julia> w = DVec(BoseFS(prod(S), 3 => 1) => 1.0)
DVec{BoseFS{1, 4, BitString{4, 1, UInt8}},Float64} with 1 entry, style = IsDeterministic{Float64}()
  fs"|0 0 1 0⟩" => 1.0

julia> dot(w, Lz, v)
0.0 + 1.0im

source

Hamiltonians interface

Behind the implementation of a particular model is a more abstract interface for defining Hamiltonians. If you want to define a new model you should make use of this interface. The most general form of a model Hamiltonian should subtype to AbstractHamiltonian and implement the relevant methods.

Rimu.Interfaces.AbstractHamiltonian — Type

AbstractHamiltonian{T}

Supertype that provides an interface for linear operators over a linear space with scalar type T that are suitable for FCIQMC (with lomc!). Indexing is done with addresses (typically not integers) from an address space that may be large (and will not need to be completely generated).

AbstractHamiltonian instances operate on vectors of type AbstractDVec from the module DictVectors and work well with addresses of type AbstractFockAddress from the module BitStringAddresses. The type works well with the external package KrylovKit.jl.

For available implementations see Hamiltonians.

Interface

Methods that need to be implemented:

Optional methods to implement:

LOStructure(::Type{typeof(lo)}): defaults to AdjointUnknown
dimension(::AbstractHamiltonian, addr): defaults to dimension of address space
allowed_address_type(h::AbstractHamiltonian): defaults to typeof(starting_address(h))
momentum(::AbstractHamiltonian): no default

Provides:

offdiagonals: iterator over reachable off-diagonal matrix elements
random_offdiagonal: function to generate random off-diagonal matrix element
*(H, v): deterministic matrix-vector multiply (allocating)
H(v): equivalent to H * v.
mul!(w, H, v): mutating matrix-vector multiply.
dot(x, H, v): compute x⋅(H*v) minimizing allocations.
H[address1, address2]: indexing with getindex() - mostly for testing purposes (slow!)
BasisSetRep: construct a basis set repesentation
sparse, Matrix: construct a (sparse) matrix representation

Geometry

Lattices in higher dimensions are defined here for HubbardRealSpace.

Rimu.Hamiltonians.LatticeGeometry — Type

abstract type LatticeGeometry{D}

A LatticeGeometry controls which sites in an AbstractFockAddress are considered to be neighbours.

Currently only supported by HubbardRealSpace.

Available implementations

Interface to implement

Base.size: return the lattice size.
neighbour_site(::LatticeGeometry, ::Int, ::Int)
num_dimensions(::LatticeGeometry)
num_neighbours(::LatticeGeometry)

source

Rimu.Hamiltonians.PeriodicBoundaries — Type

PeriodicBoundaries(size...) <: LatticeGeometry

Rectangular lattice with periodic boundary conditions of size size.

The dimension of the lattice is controlled by the number of arguments given to its constructor.

This is the default geometry used by HubbardRealSpace.

Example

julia> lattice = PeriodicBoundaries(5, 4) # 2D lattice of size 5 × 4
PeriodicBoundaries(5, 4)

julia> num_neighbours(lattice)
4

julia> neighbour_site(lattice, 1, 1)
2

julia> neighbour_site(lattice, 1, 2)
5

julia> neighbour_site(lattice, 1, 3)
6

julia> neighbour_site(lattice, 1, 4)
16

See also

source

Rimu.Hamiltonians.HardwallBoundaries — Type

HardwallBoundaries

Rectangular lattice with hard wall boundary conditions of size size. neighbour_site() will return 0 for some neighbours of boundary sites.

The dimension of the lattice is controlled by the number of arguments given to its constructor.

Example

julia> lattice = HardwallBoundaries(5) # 1D lattice of size 5
HardwallBoundaries(5)

julia> neighbour_site(lattice, 1, 1)
2

julia> neighbour_site(lattice, 1, 2)
0

See also

source

Rimu.Hamiltonians.LadderBoundaries — Type

LadderBoundaries(size...; subgeometry=PeriodicBoundaries) <: LatticeGeometry

Lattice geometry where the first dimension is of size 2 and has hardwall boundary conditions. Using this geometry is more efficient than using HardwallBoundaries with a size of 2, as it does not generate rejected neighbours.

In other dimensions, it behaves like its subgeometry, which can be any LatticeGeometry.

Example

julia> lattice = LadderBoundaries(2, 3, 4) # 3D lattice of size 2 × 3 × 4
LadderBoundaries(2, 3, 4)

julia> num_neighbours(lattice)
5

julia> neighbour_site(lattice, 1, 1)
2

julia> neighbour_site(lattice, 1, 2)
3

julia> neighbour_site(lattice, 1, 3)
5

julia> neighbour_site(lattice, 1, 4)
7

julia> neighbour_site(lattice, 1, 5)
19

See also

source

Rimu.Hamiltonians.num_neighbours — Function

num_neighbours(geom::LatticeGeometry)

Return the number of neighbours each lattice site has in this geometry.

Note that for efficiency reasons, all sites are expected to have the same number of neighbours. If some of the neighbours are invalid, this is handled by having neighbour_site return 0.

source

Rimu.Hamiltonians.num_dimensions — Function

num_dimensions(geom::LatticeGeometry)

Return the number of dimensions of the lattice in this geometry.

source

Rimu.Hamiltonians.neighbour_site — Function

neighbour_site(geom::LatticeGeometry, site, i)

Find the i-th neighbour of site in the geometry. If the move is illegal, return 0.

source

Harmonic Oscillator

Useful utilities for harmonic oscillator in Cartesian basis, see HOCartesianContactInteractions and HOCartesianEnergyConservedPerDim.

Rimu.Hamiltonians.get_all_blocks — Function

get_all_blocks(h::Union{HOCartesianContactInteractions,HOCartesianEnergyConservedPerDim}; 
    target_energy = nothing, 
    max_energy = nothing, 
    max_blocks = nothing, 
    method = :vertices,
    kwargs...) -> df

Find all distinct blocks of h. Returns a DataFrame with columns

block_id: index of block in order found
block_E0: noninteracting energy of all elements in the block
block_size: number of elements in the block
addr: first address that generates the block with e.g. BasisSetRep
indices: tuple of mode indices that allow recreation of the generating address addr; in this case use e.g. BoseFS(M; indices .=> 1) This is useful when the DataFrame is loaded from file since Arrow.jl converts custom types to NamedTuples.
t_basis: time to generate the basis for each block

Keyword arguments:

target_energy: only blocks with this noninteracting energy are found
max_energy: only blocks with noninteracting energy less than this are found
max_blocks: exit after finding this many blocks
method: Choose between :vertices and :comb for method of enumerating tuples of quantum numbers
save_to_file=nothing: if set then the DataFrame recording blocks is saved after each new block is found
additional kwargs: passed to isapprox for comparing block energies. Useful for anisotropic traps

Note: If h was constructed with option block_by_level = false then the block seeds 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 with an appropriate filter to generate the blocks.

source

Rimu.Hamiltonians.fock_to_cart — Function

fock_to_cart(addr, S; zero_index = true)

Convert a Fock state address addr to Cartesian harmonic oscillator basis indices $n_x,n_y,\ldots$. These indices are bounded by S which is a tuple of the maximum number of states in each dimension. By default the groundstate in each dimension is indexed by 0, but this can be changed by setting zero_index = false.

source

Underlying integrals for the interaction matrix elements are implemented in the following unexported functions

Rimu.Hamiltonians.four_oscillator_integral_general — Function

four_oscillator_integral_general(i, j, k, l; max_level = typemax(Int))

Integral of four one-dimensional harmonic oscillator functions,

\[ \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. Titchmarsh (1948). This is a generalisation of the closed form in Papenbrock (2002), which is is the special case where $i+j == k+l$, but is numerically unstable for large arguments. Used in HOCartesianContactInteractions and HOCartesianEnergyConservedPerDim.

source

Rimu.Hamiltonians.ho_delta_potential — Function

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,

\[ v_{ij} = \phi_{\mathbf{n}_i}(0) \phi_{\mathbf{n}_j}(0)\]

which is only non-zero for even-parity states. The ith 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.