From 5f48b35104bc636c059ed106be8110d5d05c1826 Mon Sep 17 00:00:00 2001 From: Joachim Brand Date: Fri, 13 Dec 2024 08:30:39 +1300 Subject: [PATCH] Interface tests - remove two component models (#299) # Remove legacy two-component address type and models This PR removes the address type `BoseFS2C`, which has been superseded by `CompositeFS` and model Hamiltonians that depend on it. A new package [`RimuLegacyHamiltonians.jl`](https://github.com/RimuQMC/RimuLegacyHamiltonians.jl) has been created to house this functionality. ## Breaking changes - remove `BoseFS2C`, `BoseHubbardMom1D2C`, `BoseHubbardReal1D2C` ## New features - new module `Rimu.InterfaceTests` with functions - `test_observable_interface` - `test_operator_interface` - `test_hamiltonian_interface` - `test_hamiltonian_structure` ## Documentation changes - The previous page on the `Hamiltonians` module in the "Developer Documentation" has been split into "Hamiltonians" and "Custom Hamiltonians" and moved into the "User Documentation" - various minor changes to the documentation and docstrings, including additional links. --------- Co-authored-by: mtsch --- Project.toml | 4 +- benchmark/benchmarks.jl | 10 - docs/make.jl | 6 +- docs/src/custom_hamiltonians.md | 110 ++++++ docs/src/hamiltonians.md | 80 +---- docs/src/projectormontecarlo.md | 12 +- docs/src/randomnumbers.md | 12 - docs/src/testing.md | 6 +- src/BitStringAddresses/multicomponent.jl | 47 --- src/Hamiltonians/BoseHubbardMom1D2C.jl | 144 -------- src/Hamiltonians/BoseHubbardReal1D2C.jl | 136 -------- src/Hamiltonians/DensityMatrixDiagonal.jl | 12 - src/Hamiltonians/G2MomCorrelator.jl | 77 +++++ src/Hamiltonians/Hamiltonians.jl | 12 +- src/Hamiltonians/HubbardRealSpace.jl | 4 +- src/Hamiltonians/Momentum.jl | 4 - src/Hamiltonians/TRSymmetry.jl | 7 - src/Hamiltonians/abstract.jl | 6 - src/Hamiltonians/correlation_functions.jl | 160 --------- src/Hamiltonians/geometry.jl | 24 +- src/Hamiltonians/ho-cart-tools.jl | 81 ++--- src/InterfaceTests.jl | 285 ++++++++++++++++ src/Rimu.jl | 2 + src/strategies_and_params/poststepstrategy.jl | 2 +- test/BitStringAddresses.jl | 18 - test/Hamiltonians.jl | 313 +----------------- test/allocation_tests.jl | 8 - test/lomc.jl | 6 +- test/runtests.jl | 45 --- 29 files changed, 578 insertions(+), 1055 deletions(-) create mode 100644 docs/src/custom_hamiltonians.md delete mode 100644 docs/src/randomnumbers.md delete mode 100644 src/Hamiltonians/BoseHubbardMom1D2C.jl delete mode 100644 src/Hamiltonians/BoseHubbardReal1D2C.jl create mode 100644 src/Hamiltonians/G2MomCorrelator.jl create mode 100644 src/InterfaceTests.jl diff --git a/Project.toml b/Project.toml index aac26b5d6..d0271dbe3 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Rimu" uuid = "c53c40cc-bd84-11e9-2cf4-a9fde2b9386e" authors = ["Joachim Brand "] -version = "0.13.2-dev" +version = "0.14.0" [deps] Arrow = "69666777-d1a9-59fb-9406-91d4454c9d45" @@ -39,6 +39,7 @@ StrLiterals = "68059f60-971f-57ff-a2d0-18e7de9ccc84" TOML = "fa267f1f-6049-4f14-aa54-33bafae1ed76" Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c" TerminalLoggers = "5d786b92-1e48-4d6f-9151-6b4477ca9bed" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" TupleTools = "9d95972d-f1c8-5527-a6e0-b4b365fa01f6" VectorInterface = "409d34a3-91d5-4945-b6ec-7529ddf182d8" @@ -94,7 +95,6 @@ julia = "1.9" Arpack = "7d9fca2a-8960-54d3-9f78-7d1dccf2cb97" IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153" KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77" -Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] test = ["Test"] diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl index f5d0e1947..4631b2832 100644 --- a/benchmark/benchmarks.jl +++ b/benchmark/benchmarks.jl @@ -63,16 +63,6 @@ const SUITE = @benchmarkset "Rimu" begin lomc!(ham, dv; s_strat, post_step, dτ=1e-4, laststep=8000) end seconds=150 - @case "(4+1, 11) 2C Mom space with G2Correlators" begin - addr = BoseFS2C(ntuple(i -> ifelse(i == 5, 4, 0), 11), ntuple(==(5), 11)) - ham = BoseHubbardMom1D2C(addr, v=0.1) - dv = PDVec(addr => 1.0f0; style=IsDynamicSemistochastic{Float32}()) - s_strat = DoubleLogUpdate(target_walkers=10_000) - replica_strategy = AllOverlaps(2; operator = ntuple(i -> G2Correlator(i - 1), 7)) - - lomc!(ham, dv; s_strat, replica_strategy, laststep=2000) - end seconds=150 - @case "(50, 50) Real space" begin addr = near_uniform(BoseFS{50,50}) ham = HubbardReal1D(addr, u=6.0) diff --git a/docs/make.jl b/docs/make.jl index f64e1adce..f06912a93 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -43,7 +43,7 @@ for fn in EXAMPLES_FILES end makedocs(; - modules=[Rimu,Rimu.RimuIO], + modules=[Rimu,Rimu.RimuIO,Rimu.InterfaceTests], format=Documenter.HTML( prettyurls = false, size_threshold=700_000, # 700 kB @@ -55,17 +55,17 @@ makedocs(; "User documentation" => [ "Exact Diagonalization" => "exactdiagonalization.md", "Projector Monte Carlo" => "projectormontecarlo.md", + "Hamiltonians" => "hamiltonians.md", "StatsTools" => "statstools.md", "Using MPI" => "mpi.md", + "Custom Hamiltonians" => "custom_hamiltonians.md", ], "Developer documentation" => [ "Interfaces" => "interfaces.md", - "Hamiltonians" => "hamiltonians.md", "Dict vectors" => "dictvectors.md", "BitString addresses" => "addresses.md", "Stochastic styles" => "stochasticstyles.md", "I/O" => "rimuio.md", - "Random numbers" => "randomnumbers.md", "Documentation generation" => "documentation.md", "Code testing" => "testing.md", ], diff --git a/docs/src/custom_hamiltonians.md b/docs/src/custom_hamiltonians.md new file mode 100644 index 000000000..825e4c5c1 --- /dev/null +++ b/docs/src/custom_hamiltonians.md @@ -0,0 +1,110 @@ +# Advanced operator usage and custom Hamiltonians + +`Rimu` can be used to work with custom Hamiltonians and observables that are user-defined and + not part of the `Rimu.jl` package. To make this possible and reliable, `Rimu` exposes a number + of interfaces and provides helper functions to test compliance with the interfaces through the + submodule [`Rimu.InterfaceTests`](@ref), see [Interface tests](@ref). This section covers the + relevant interfaces, the interface functions as well as potentially useful helper functions. + + In order to define custom Hamiltonians or observables it is useful to know how the operator + type hierarchy works in `Rimu`. For an example of how to implement custom Hamiltonians that + are not part of the `Rimu.jl` package, see + [`RimuLegacyHamiltonians.jl`](https://github.com/RimuQMC/RimuLegacyHamiltonians.jl). + +## Operator type hierarchy + +`Rimu` offers a hierarchy of abstract types that define interfaces with different requirements +for operators: +```julia +AbstractHamiltonian <: AbstractOperator <: AbstractObservable +``` +The different abstract types have different requirements and are meant to be used for different purposes. +- [`AbstractHamiltonian`](@ref)s are fully featured models that define a Hilbert space and a linear operator over a scalar field. They can be passed as a Hamiltonian into [`ProjectorMonteCarloProblem`](@ref) or [`ExactDiagonalizationProblem`](@ref). +- [`AbstractOperator`](@ref) and [`AbstractObservable`](@ref) are supertypes of [`AbstractHamiltonian`](@ref) with less stringent conditions. They are useful for defining observables that can be used in a three-way `dot` product, or passed as observables into a [`ReplicaStrategy`](@ref) that can be inserted with the keyword `replica_strategy` into a [`ProjectorMonteCarloProblem`](@ref). + +## 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. A new +model Hamiltonian should subtype to `AbstractHamiltonian` and implement the relevant methods. + +```@docs +AbstractHamiltonian +offdiagonals +diagonal_element +starting_address +``` + +The following functions may be implemented instead of [`offdiagonals`](@ref). + +```@docs +num_offdiagonals +get_offdiagonal +``` + +The following functions come with default implementations, but may be customized. + +```@docs +random_offdiagonal +Hamiltonians.LOStructure +dimension +has_adjoint +allows_address_type +Base.eltype +VectorInterface.scalartype +mul! +``` + +This interface relies on unexported functionality, including +```@docs +Hamiltonians.adjoint +Hamiltonians.dot +Hamiltonians.AbstractOffdiagonals +Hamiltonians.Offdiagonals +Hamiltonians.check_address_type +Hamiltonians.number_conserving_dimension +Hamiltonians.number_conserving_bose_dimension +Hamiltonians.number_conserving_fermi_dimension +``` + +## Operator and observable interface + +```@docs +AbstractObservable +AbstractOperator +``` + +## Interface tests +Helper functions that can be used for testing the various interfaces are provided in the +(unexported) submodule `Rimu.InterfaceTests`. + +```@docs +Rimu.InterfaceTests +``` + +### Testing functions +```@docs +Rimu.InterfaceTests.test_hamiltonian_interface +Rimu.InterfaceTests.test_hamiltonian_structure +Rimu.InterfaceTests.test_observable_interface +Rimu.InterfaceTests.test_operator_interface +``` + +## Utilities for harmonic oscillator models +Useful utilities for harmonic oscillator in Cartesian basis, see [`HOCartesianContactInteractions`](@ref) +and [`HOCartesianEnergyConservedPerDim`](@ref). +```@docs +get_all_blocks +fock_to_cart +``` +Underlying integrals for the interaction matrix elements are implemented in the following unexported functions +```@docs +Hamiltonians.four_oscillator_integral_general +Hamiltonians.ho_delta_potential +Hamiltonians.log_abs_oscillator_zero +``` + +## Index +```@index +Pages = ["custom_hamiltonians.md"] +``` diff --git a/docs/src/hamiltonians.md b/docs/src/hamiltonians.md index 94d3855d5..d1d9eb081 100644 --- a/docs/src/hamiltonians.md +++ b/docs/src/hamiltonians.md @@ -11,44 +11,39 @@ The Hamiltonians can be used for projector quantum Monte Carlo with [`ProjectorM Hamiltonians ``` - -## 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 +## Real space Hubbard models ```@docs HubbardReal1D -BoseHubbardReal1D2C HubbardReal1DEP HubbardRealSpace ExtendedHubbardReal1D ``` -### Momentum space Hubbard models +## Momentum space Hubbard models ```@docs HubbardMom1D -BoseHubbardMom1D2C HubbardMom1DEP ExtendedHubbardMom1D ``` -### Harmonic oscillator models +## Harmonic oscillator models ```@docs HOCartesianContactInteractions HOCartesianEnergyConservedPerDim HOCartesianCentralImpurity ``` -### Other +## Other model Hamiltonians ```@docs MatrixHamiltonian Transcorrelated1D FroehlichPolaron ``` -### Convenience functions +## Convenience functions ```@docs rayleigh_quotient momentum @@ -81,8 +76,6 @@ AbstractHamiltonian{T} <: AbstractOperator{T} <: AbstractObservable{T} ``` ```@docs -AbstractObservable -AbstractOperator ParticleNumberOperator G2RealCorrelator G2RealSpace @@ -97,55 +90,10 @@ Momentum AxialAngularMomentumHO ``` -## 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. - -```@docs -AbstractHamiltonian -offdiagonals -diagonal_element -starting_address -``` - -The following functions may be implemented instead of [`offdiagonals`](@ref). - -```@docs -num_offdiagonals -get_offdiagonal -``` - -The following functions come with default implementations, but may be customized. - -```@docs -random_offdiagonal -Hamiltonians.LOStructure -dimension -has_adjoint -allows_address_type -Base.eltype -VectorInterface.scalartype -mul! -``` - -This interface relies on unexported functionality, including -```@docs -Hamiltonians.adjoint -Hamiltonians.dot -Hamiltonians.AbstractOffdiagonals -Hamiltonians.Offdiagonals -Hamiltonians.check_address_type -Hamiltonians.number_conserving_dimension -Hamiltonians.number_conserving_bose_dimension -Hamiltonians.number_conserving_fermi_dimension -``` - ## Geometry -Lattices in higher dimensions are defined here for [`HubbardRealSpace`](@ref) and [`G2RealSpace`](@ref). +Lattices in higher dimensions are defined here and can be passed with the keyword argument +`geometry` to [`HubbardRealSpace`](@ref) and [`G2RealSpace`](@ref). ```@docs CubicGrid @@ -157,20 +105,6 @@ HardwallBoundaries LadderBoundaries ``` -## Harmonic Oscillator -Useful utilities for harmonic oscillator in Cartesian basis, see [`HOCartesianContactInteractions`](@ref) -and [`HOCartesianEnergyConservedPerDim`](@ref). -```@docs -get_all_blocks -fock_to_cart -``` -Underlying integrals for the interaction matrix elements are implemented in the following unexported functions -```@docs -Hamiltonians.four_oscillator_integral_general -Hamiltonians.ho_delta_potential -Hamiltonians.log_abs_oscillator_zero -``` - ## Index ```@index Pages = ["hamiltonians.md"] diff --git a/docs/src/projectormontecarlo.md b/docs/src/projectormontecarlo.md index db26f5d8b..64ab9da7a 100644 --- a/docs/src/projectormontecarlo.md +++ b/docs/src/projectormontecarlo.md @@ -7,9 +7,9 @@ Full Configuration Interaction Quantum Monte Carlo (FCIQMC). ## `ProjectorMonteCarloProblem` -To run a projector Monte Carlo simulation you set up a problem with `ProjectorMonteCarloProblem` -and solve it with `solve`. Alternatively you can initialize a `PMCSimulation` struct, `step!` -through time steps, and `solve!` it to completion. +To run a projector Monte Carlo simulation you set up a problem with [`ProjectorMonteCarloProblem`](@ref) +and solve it with [`solve`](@ref). Alternatively you can [`init`](@ref) it with to obtain a [`PMCSimulation`](@ref Rimu.PMCSimulation) struct, [`step!`](@ref) +through time steps, and [`solve!`](@ref) it to completion. ```@docs; canonical=false ProjectorMonteCarloProblem @@ -22,16 +22,16 @@ step! After `solve` or `solve!` have been called the returned `PMCSimulation` contains the results of the projector Monte Carlo calculation. -### `PMCSimulation` and report as a `DataFrame` +## `PMCSimulation` and report as a `DataFrame` ```@docs; canonical=false Rimu.PMCSimulation ``` -The `DataFrame` returned from `DataFrame(::PMCSimulation)` contains the time series data from +The [`DataFrame`](https://dataframes.juliadata.org/stable/) returned from `DataFrame(::PMCSimulation)` contains the time series data from the projector Monte Carlo simulation that is of primary interest for analysis. Depending on the `reporting_strategy` and other options passed as keyword arguments to -`ProjectorMonteCarloProblem` it can have different numbers of rows and columns. The rows +[`ProjectorMonteCarloProblem`](@ref) it can have different numbers of rows and columns. The rows correspond to the reported time steps (Monte Carlo steps). There is at least one column with the name `:step`. Further columns are usually present with additional data reported from the simulation. For the default option `algorithm = FCIQMC(; shift_strategy, time_step_strategy)` with a single diff --git a/docs/src/randomnumbers.md b/docs/src/randomnumbers.md deleted file mode 100644 index fb6d29207..000000000 --- a/docs/src/randomnumbers.md +++ /dev/null @@ -1,12 +0,0 @@ -# Random numbers in Rimu - -Rimu uses Julia's built-in random number generator, which currently defaults to -[Xoshiro256++](https://docs.julialang.org/en/v1/stdlib/Random/#Random.Xoshiro). - -## Reproducibility - -If you want FCIQMC runs to be reproducible, make sure to seed the RNG with -[Random.seed!](https://docs.julialang.org/en/v1/stdlib/Random/#Random.seed!). - -MPI-distributed runs can also be made reproducible by seeding the RNG with -[`mpi_seed!`](@ref). diff --git a/docs/src/testing.md b/docs/src/testing.md index 54ba295eb..761a286cf 100644 --- a/docs/src/testing.md +++ b/docs/src/testing.md @@ -1,6 +1,6 @@ # Code testing -The script `runtest.jl` in the `test/` folder contains tests of the code. To run the test simply run the script from the Julia REPL or run +The script `runtest.jl` in the `test/` folder contains tests of the code in `Rimu`. To run the test simply run the script from the Julia REPL or run ``` Rimu$ julia test/runtest.jl ``` @@ -12,3 +12,7 @@ More tests should be added over time to test core functionality of the code. To GitHub Actions are set up to run the test script automatically on the GitHub cloud server every time a new commit to the master branch is pushed to the server. The setup for this to happen is configured in the file `actions.yml` in the `Rimu/.github/workflows` folder. + +## Testing of custom types for use with `Rimu` + +The module `Rimu.InterfaceTests` contains a number of functions to test the interfaces of the [`AbstractHamiltonian`](@ref) type hierarchy. See [Interface tests](@ref) in the section [Advanced operator usage and custom Hamiltonians](@ref). diff --git a/src/BitStringAddresses/multicomponent.jl b/src/BitStringAddresses/multicomponent.jl index d67b2796e..e8cf23cac 100644 --- a/src/BitStringAddresses/multicomponent.jl +++ b/src/BitStringAddresses/multicomponent.jl @@ -1,47 +1,3 @@ -""" - BoseFS2C{NA,NB,M,AA,AB} <: AbstractFockAddress - BoseFS2C(onr_a, onr_b) - -Address type that constructed with two [`BoseFS{N,M,S}`](@ref). It represents a -Fock state with two components, e.g. two different species of bosons with particle -number `NA` from species S and particle number `NB` from species B. The number of -modes `M` is expected to be the same for both components. -""" -struct BoseFS2C{NA,NB,M,SA,SB,N} <: AbstractFockAddress{N,M} - bsa::BoseFS{NA,M,SA} - bsb::BoseFS{NB,M,SB} - -end - -function BoseFS2C(bsa::BoseFS{NA,M,SA}, bsb::BoseFS{NB,M,SB}) where {NA,NB,M,SA,SB} - N = NA + NB - return BoseFS2C{NA,NB,M,SA,SB,N}(bsa, bsb) -end -BoseFS2C(onr_a::Tuple, onr_b::Tuple) = BoseFS2C(BoseFS(onr_a),BoseFS(onr_b)) - -function print_address(io::IO, b::BoseFS2C; compact) - if compact - print_address(io, b.bsa; compact) - print(io, " ⊗ ") - print_address(io, b.bsb; compact) - else - print(io, "BoseFS2C(", b.bsa, ", ", b.bsb, ")") - end -end - -num_components(::Type{<:BoseFS2C}) = 2 -Base.isless(a::T, b::T) where {T<:BoseFS2C} = isless((a.bsa, a.bsb), (b.bsa, b.bsb)) -onr(b2::BoseFS2C) = (onr(b2.bsa), onr(b2.bsb)) - - -function near_uniform(::Type{<:BoseFS2C{NA,NB,M}}) where {NA,NB,M} - return BoseFS2C(near_uniform(BoseFS{NA,M}), near_uniform(BoseFS{NB,M})) -end - -function time_reverse(c::BoseFS2C{NA,NA,M,S,S,N}) where {NA,M,S,N} - return BoseFS2C{NA,NA,M,S,S,N}(c.bsb, c.bsa) -end - """ CompositeFS(addresses::SingleComponentFockAddress...) <: AbstractFockAddress @@ -190,6 +146,3 @@ function print_address(io::IO, f::FermiFS2C; compact=false) invoke(print_address, Tuple{typeof(io),CompositeFS}, io, f) end end - -BoseFS2C(fs::CompositeFS{2}) = BoseFS2C(fs.components...) -CompositeFS(fs::BoseFS2C) = CompositeFS(fs.bsa, fs.bsb) diff --git a/src/Hamiltonians/BoseHubbardMom1D2C.jl b/src/Hamiltonians/BoseHubbardMom1D2C.jl deleted file mode 100644 index 76c12c91e..000000000 --- a/src/Hamiltonians/BoseHubbardMom1D2C.jl +++ /dev/null @@ -1,144 +0,0 @@ -""" - 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. - -```math -\\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 - -* `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. -* `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`](@ref). - -# See also - -* [`BoseFS2C`](@ref) -* [`BoseHubbardReal1D2C`](@ref) - -""" -struct BoseHubbardMom1D2C{T,HA,HB,V} <: TwoComponentHamiltonian{T} - ha::HA - hb::HB -end - -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 - -function Base.show(io::IO, h::BoseHubbardMom1D2C) - addr = starting_address(h) - ua = h.ha.u - ub = h.hb.u - ta = h.ha.t - tb = h.hb.t - v = h.v - print(io, "BoseHubbardMom1D2C($addr; ua=$ua, ub=$ub, ta=$ta, tb=$tb, v=$v)") -end - -dimension(::BoseHubbardMom1D2C, address) = number_conserving_dimension(address) - -function starting_address(h::BoseHubbardMom1D2C) - return BoseFS2C(starting_address(h.ha), starting_address(h.hb)) -end - -LOStructure(::Type{<:BoseHubbardMom1D2C{<:Real}}) = IsHermitian() - -Base.getproperty(h::BoseHubbardMom1D2C, s::Symbol) = getproperty(h, Val(s)) -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, 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, 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 - for k in 1:M - interaction2c += onrep_a[k]*onrep_b[p] # b†_p b_p a†_k a_k - end - end - return ( - diagonal_element(ham.ha, address.bsa) + - diagonal_element(ham.hb, address.bsb) + - ham.v/M*interaction2c - ) -end - -function get_offdiagonal(ham::BoseHubbardMom1D2C, address::BoseFS2C, chosen) - return offdiagonals(ham, address)[chosen] -end - -""" - OffdiagonalsBoseMom1D2C - -Specialized [`AbstractOffdiagonals`](@ref) that keep track of number of off-diagonals and -number of occupied sites in both components of the address. -""" -struct OffdiagonalsBoseMom1D2C{ - A<:BoseFS2C,T,H<:TwoComponentHamiltonian{T},O1,O2,M1,M2 -} <: AbstractOffdiagonals{A,T} - hamiltonian::H - address::A - length::Int - offdiags_a::O1 - offdiags_b::O2 - map_a::M1 - map_b::M2 -end - -function offdiagonals(h::BoseHubbardMom1D2C{T}, a::BoseFS2C) where {T} - offdiags_a = offdiagonals(h.ha, a.bsa) - offdiags_b = offdiagonals(h.hb, a.bsb) - map_a = OccupiedModeMap(a.bsa) - map_b = OccupiedModeMap(a.bsb) - occ_a = length(map_a) - occ_b = length(map_b) - len = length(offdiags_a) + length(offdiags_b) + occ_a * (num_modes(a) - 1) * occ_b - - return OffdiagonalsBoseMom1D2C(h, a, len, offdiags_a, offdiags_b, map_a, map_b) -end - -function Base.getindex(s::OffdiagonalsBoseMom1D2C{A,T}, i)::Tuple{A,T} where {A,T} - @boundscheck 1 ≤ i ≤ s.length || throw(BoundsError(s, i)) - num_hops_a = length(s.offdiags_a) - num_hops_b = length(s.offdiags_b) - if i ≤ num_hops_a - new_a, matrix_element = s.offdiags_a[i] - new_add = A(new_a, s.address.bsb) - elseif i ≤ num_hops_a + num_hops_b - i -= num_hops_a - new_b, matrix_element = s.offdiags_b[i] - new_add = A(s.address.bsa, new_b) - else - i -= num_hops_a + num_hops_b - new_a, new_b, val = momentum_transfer_excitation( - s.address.bsa, s.address.bsb, i, s.map_a, s.map_b - ) - new_add = A(new_a, new_b) - matrix_element = s.hamiltonian.v/num_modes(new_add) * val - end - return new_add, matrix_element -end - -Base.size(s::OffdiagonalsBoseMom1D2C) = (s.length,) diff --git a/src/Hamiltonians/BoseHubbardReal1D2C.jl b/src/Hamiltonians/BoseHubbardReal1D2C.jl deleted file mode 100644 index 74411e6d4..000000000 --- a/src/Hamiltonians/BoseHubbardReal1D2C.jl +++ /dev/null @@ -1,136 +0,0 @@ -""" - 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. - -```math -\\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 - -* [`HubbardReal1D`](@ref) -* [`BoseHubbardMom1D2C`](@ref) - -""" -struct BoseHubbardReal1D2C{T,HA,HB,V} <: TwoComponentHamiltonian{T} - ha::HA - hb::HB -end - -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 - -function Base.show(io::IO, h::BoseHubbardReal1D2C) - addr = starting_address(h) - ua = h.ha.u - ub = h.hb.u - ta = h.ha.t - tb = h.hb.t - v = h.v - print(io, "BoseHubbardReal1D2C($addr; ua=$ua, ub=$ub, ta=$ta, tb=$tb, v=$v)") -end - -dimension(::BoseHubbardReal1D2C, address) = number_conserving_dimension(address) - -function starting_address(h::BoseHubbardReal1D2C) - return BoseFS2C(starting_address(h.ha), starting_address(h.hb)) -end - -LOStructure(::Type{<:BoseHubbardReal1D2C{<:Real}}) = IsHermitian() - -Base.getproperty(h::BoseHubbardReal1D2C, s::Symbol) = getproperty(h, Val(s)) -Base.getproperty(h::BoseHubbardReal1D2C, ::Val{:ha}) = getfield(h, :ha) -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(::BoseHubbardReal1D2C, address) - return 2*(num_occupied_modes(address.bsa) + num_occupied_modes(address.bsb)) -end - -""" - bose_hubbard_2c_interaction(::BoseFS2C) - -Compute the interaction between the two components. -""" -function bose_hubbard_2c_interaction(address::BoseFS2C) - c1 = onr(address.bsa) - c2 = onr(address.bsb) - interaction = zero(eltype(c1)) - for site = 1:length(c1) - if !iszero(c2[site]) - interaction += c2[site] * c1[site] - end - end - return interaction -end - -function diagonal_element(ham::BoseHubbardReal1D2C, address::BoseFS2C) - return ( - diagonal_element(ham.ha, address.bsa) + - diagonal_element(ham.hb, address.bsb) + - ham.v * bose_hubbard_2c_interaction(address) - ) -end - -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(address.bsa, chosen) - elem = - ham.ha.t * onproduct - return BoseFS2C(naddress_from_bsa,address.bsb), elem - else - chosen -= nhops_a - naddress_from_bsb, onproduct = hopnextneighbour(address.bsb, chosen) - elem = -ham.hb.t * onproduct - return BoseFS2C(address.bsa,naddress_from_bsb), elem - end - # return new address and matrix element -end - -struct OffdiagonalsBoseReal1D2C{ - A<:BoseFS2C,T,H<:TwoComponentHamiltonian{T} -} <: AbstractOffdiagonals{A,T} - hamiltonian::H - address::A - length::Int - num_hops_a::Int -end - -function offdiagonals(h::BoseHubbardReal1D2C, a::BoseFS2C) - hops_a = num_offdiagonals(h.ha, a.bsa) - hops_b = num_offdiagonals(h.hb, a.bsb) - length = hops_a + hops_b - - return OffdiagonalsBoseReal1D2C(h, a, length, hops_a) -end - -function Base.getindex(s::OffdiagonalsBoseReal1D2C{A,T}, i)::Tuple{A,T} where {A,T} - @boundscheck 1 ≤ i ≤ s.length || throw(BoundsError(s, i)) - if i ≤ s.num_hops_a - new_a, matrix_element = get_offdiagonal(s.hamiltonian.ha, s.address.bsa, i) - new_add = A(new_a, s.address.bsb) - else - i -= s.num_hops_a - new_b, matrix_element = get_offdiagonal(s.hamiltonian.hb, s.address.bsb, i) - new_add = A(s.address.bsa, new_b) - end - return new_add, matrix_element -end - -Base.size(s::OffdiagonalsBoseReal1D2C) = (s.length,) diff --git a/src/Hamiltonians/DensityMatrixDiagonal.jl b/src/Hamiltonians/DensityMatrixDiagonal.jl index b4f84e08e..e63de7b29 100644 --- a/src/Hamiltonians/DensityMatrixDiagonal.jl +++ b/src/Hamiltonians/DensityMatrixDiagonal.jl @@ -40,17 +40,5 @@ function diagonal_element(dmd::DensityMatrixDiagonal{C}, add::CompositeFS) where return float(find_mode(comp, dmd.mode).occnum) end -function diagonal_element(dmd::DensityMatrixDiagonal{0}, add::BoseFS2C) - return float(find_mode(add.bsa, dmd.mode).occnum + find_mode(add.bsb, dmd.mode).occnum) -end -function diagonal_element(dmd::DensityMatrixDiagonal{1}, add::BoseFS2C) - comp = add.bsa - return float(find_mode(comp, dmd.mode).occnum) -end -function diagonal_element(dmd::DensityMatrixDiagonal{2}, add::BoseFS2C) - comp = add.bsb - return float(find_mode(comp, dmd.mode).occnum) -end - num_offdiagonals(dmd::DensityMatrixDiagonal, _) = 0 LOStructure(::Type{<:DensityMatrixDiagonal}) = IsDiagonal() diff --git a/src/Hamiltonians/G2MomCorrelator.jl b/src/Hamiltonians/G2MomCorrelator.jl new file mode 100644 index 000000000..08534d389 --- /dev/null +++ b/src/Hamiltonians/G2MomCorrelator.jl @@ -0,0 +1,77 @@ + +import Rimu.Hamiltonians: num_offdiagonals, diagonal_element, get_offdiagonal + +""" + G2MomCorrelator(d::Int) <: AbstractOperator{ComplexF64} + +Two-body correlation operator representing the density-density +correlation at distance `d`. It returns a `Complex` value. + +Correlation within a single component: +```math +\\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 +```math +\\frac{1}{M}\\sum_{k,p=1}^M a^†_{k} b^†_{p} b_p a_k . +``` + +# Arguments +- `d::Integer`: the distance between two particles. + + +# See also + +* [`Rimu.G2RealCorrelator`](@ref) +* [`Rimu.G2RealSpace`](@ref) +* [`Rimu.AbstractOperator`](@ref) +* [`Rimu.AllOverlaps`](@ref) +""" +struct G2MomCorrelator{C} <: AbstractOperator{ComplexF64} + d::Int +end +# The type parameter `C` is not used here, but may be used for future extensions. +# It is kept here for consistency with `RimuLegacyHamiltonians.jl`. +function G2MomCorrelator(d::Int) + return G2MomCorrelator{3}(d) +end + +function Rimu.Interfaces.allows_address_type(g2m::G2MomCorrelator, ::Type{A}) where {A} + return num_modes(A) > g2m.d && A <: SingleComponentFockAddress +end + +function Base.show(io::IO, g::G2MomCorrelator{3}) + # 3 is the default value for the type parameter + print(io, "G2MomCorrelator($(g.d))") +end + +function num_offdiagonals(g::G2MomCorrelator, addr::SingleComponentFockAddress) + m = num_modes(addr) + singlies, doublies = num_singly_doubly_occupied_sites(addr) + return singlies * (singlies - 1) * (m - 2) + doublies * (m - 1) +end + +function diagonal_element(g::G2MomCorrelator, addr::SingleComponentFockAddress) + M = num_modes(addr) + onrep = onr(addr) + gd = 0 + for p in 1:M + iszero(onrep[p]) && continue + for k in 1:M + gd += onrep[k] * onrep[p] # a†_p a_p a†_k a_k + end + end + return ComplexF64(gd / M) +end + +function get_offdiagonal( + g::G2MomCorrelator, + addr::A, + chosen, +)::Tuple{A,ComplexF64} where {A<:SingleComponentFockAddress} + M = num_modes(addr) + new_add, gamma, Δp = momentum_transfer_excitation(addr, chosen, OccupiedModeMap(addr)) + gd = exp(-im * g.d * Δp * 2π / M) * gamma + return new_add, ComplexF64(gd / M) +end diff --git a/src/Hamiltonians/Hamiltonians.jl b/src/Hamiltonians/Hamiltonians.jl index 493eb8e5b..05908b177 100644 --- a/src/Hamiltonians/Hamiltonians.jl +++ b/src/Hamiltonians/Hamiltonians.jl @@ -6,14 +6,12 @@ Hamiltonians. Real space Hubbard models - [`HubbardReal1D`](@ref) - - [`BoseHubbardReal1D2C`](@ref) - [`HubbardReal1DEP`](@ref) - [`HubbardRealSpace`](@ref) - [`ExtendedHubbardReal1D`](@ref) Momentum space Hubbard models - [`HubbardMom1D`](@ref) -- [`BoseHubbardMom1D2C`](@ref) - [`HubbardMom1DEP`](@ref) Harmonic oscillator models @@ -36,8 +34,8 @@ Other ## [Observables](#Observables) - [`ParticleNumberOperator`](@ref) - [`G2RealCorrelator`](@ref) -- [`G2RealSpace`](@ref) - [`G2MomCorrelator`](@ref) +- [`G2RealSpace`](@ref) - [`DensityMatrixDiagonal`](@ref) - [`SingleParticleExcitation`](@ref) - [`TwoParticleExcitation`](@ref) @@ -75,7 +73,6 @@ export MatrixHamiltonian export HubbardReal1D, HubbardMom1D, ExtendedHubbardReal1D, ExtendedHubbardMom1D, HubbardRealSpace export HubbardReal1DEP, shift_lattice, shift_lattice_inv export HubbardMom1DEP -export BoseHubbardMom1D2C, BoseHubbardReal1D2C export GutzwillerSampling, GuidingVectorSampling export ParitySymmetry export TimeReversalSymmetry @@ -85,9 +82,9 @@ export hubbard_dispersion, continuum_dispersion export FroehlichPolaron export ParticleNumberOperator -export G2MomCorrelator, G2RealCorrelator, G2RealSpace, SuperfluidCorrelator, DensityMatrixDiagonal, Momentum +export G2RealCorrelator, G2RealSpace, SuperfluidCorrelator, DensityMatrixDiagonal, Momentum export SingleParticleExcitation, TwoParticleExcitation, ReducedDensityMatrix -export StringCorrelator +export StringCorrelator, G2MomCorrelator export CubicGrid, PeriodicBoundaries, HardwallBoundaries, LadderBoundaries @@ -118,8 +115,6 @@ include("HubbardMom1DEP.jl") include("HubbardRealSpace.jl") include("ExtendedHubbardReal1D.jl") -include("BoseHubbardReal1D2C.jl") -include("BoseHubbardMom1D2C.jl") include("FroehlichPolaron.jl") include("GutzwillerSampling.jl") @@ -130,6 +125,7 @@ include("Stoquastic.jl") include("Transcorrelated1D.jl") include("correlation_functions.jl") +include("G2MomCorrelator.jl") include("DensityMatrixDiagonal.jl") include("reduced_density_matrix.jl") include("Momentum.jl") diff --git a/src/Hamiltonians/HubbardRealSpace.jl b/src/Hamiltonians/HubbardRealSpace.jl index dd1a75fb8..533204da2 100644 --- a/src/Hamiltonians/HubbardRealSpace.jl +++ b/src/Hamiltonians/HubbardRealSpace.jl @@ -201,9 +201,9 @@ function HubbardRealSpace( throw(ArgumentError("`t` must be a vector of length $C")) elseif length(v) ≠ C * D throw(ArgumentError("`v` must be a $C × $D matrix")) - elseif address isa BoseFS2C + elseif !(address isa SingleComponentFockAddress || address isa CompositeFS) throw(ArgumentError( - "`BoseFS2C` is not supported for this Hamiltonian, use `CompositeFS`" + "unsupported address type detected use `CompositeFS` or `<: SingleComponentFockAddress`" )) end warn_fermi_interaction(address, u) diff --git a/src/Hamiltonians/Momentum.jl b/src/Hamiltonians/Momentum.jl index f507fe71c..464bb8a50 100644 --- a/src/Hamiltonians/Momentum.jl +++ b/src/Hamiltonians/Momentum.jl @@ -53,9 +53,5 @@ end diagonal_element(m::Momentum{0}, a::SingleComponentFockAddress) = _momentum(a, m.fold) diagonal_element(m::Momentum{1}, a::SingleComponentFockAddress) = _momentum(a, m.fold) -diagonal_element(m::Momentum{0}, a::BoseFS2C) = _momentum((a.bsa, a.bsb), m.fold) -diagonal_element(m::Momentum{1}, a::BoseFS2C) = _momentum(a.bsa, m.fold) -diagonal_element(m::Momentum{2}, a::BoseFS2C) = _momentum(a.bsb, m.fold) - diagonal_element(m::Momentum{0}, a::CompositeFS) = _momentum(a.components, m.fold) diagonal_element(m::Momentum{N}, a::CompositeFS) where {N} = _momentum(a.components[N], m.fold) diff --git a/src/Hamiltonians/TRSymmetry.jl b/src/Hamiltonians/TRSymmetry.jl index 66e3a772a..d191d953d 100644 --- a/src/Hamiltonians/TRSymmetry.jl +++ b/src/Hamiltonians/TRSymmetry.jl @@ -57,12 +57,6 @@ function check_tr_address(addr::CompositeFS) throw(ArgumentError("Two component address with equal particle numbers and component types required for `TimeReversalSymmetry`.")) end end -function check_tr_address(addr::BoseFS2C{NA,NB,M,SA,SB,N}) where {NA,NB,M,SA,SB,N} - if NA ≠ NB || SA ≠ SB - throw(ArgumentError("Two component address with equal particle numbers required for `TimeReversalSymmetry`.")) - end -end - function Base.show(io::IO, h::TimeReversalSymmetry) print(io, "TimeReversalSymmetry(", h.hamiltonian, ", even=", h.even, ")") @@ -123,4 +117,3 @@ end function diagonal_element(h::TimeReversalSymmetry, add) return diagonal_element(h.hamiltonian, add) end - diff --git a/src/Hamiltonians/abstract.jl b/src/Hamiltonians/abstract.jl index 7b74f54bd..f2b7f9197 100644 --- a/src/Hamiltonians/abstract.jl +++ b/src/Hamiltonians/abstract.jl @@ -80,9 +80,6 @@ end function dimension(::Type{<:FermiFS{N,M}}) where {N,M} return number_conserving_fermi_dimension(N, M) end -function dimension(::Type{<:BoseFS2C{NA,NB,M}}) where {NA,NB,M} - return dimension(BoseFS{NA,M}) * dimension(BoseFS{NB,M}) -end function dimension(::Type{<:CompositeFS{<:Any,<:Any,<:Any,T}}) where {T} return prod(dimension, T.parameters) # This relies on an implementation detail of the Tuple type and may break in future @@ -140,9 +137,6 @@ function number_conserving_dimension(address::FermiFS) n = num_particles(address) return number_conserving_fermi_dimension(n, m) end -function number_conserving_dimension(addr::BoseFS2C) - return number_conserving_dimension(addr.bsa) * number_conserving_dimension(addr.bsb) -end function number_conserving_dimension(address::CompositeFS) return prod(number_conserving_dimension, address.components) end diff --git a/src/Hamiltonians/correlation_functions.jl b/src/Hamiltonians/correlation_functions.jl index dd8d08b1f..65bf247cf 100644 --- a/src/Hamiltonians/correlation_functions.jl +++ b/src/Hamiltonians/correlation_functions.jl @@ -235,166 +235,6 @@ function diagonal_element(g2::G2RealSpace{0,0}, addr::CompositeFS) return _g2_diagonal_element(g2, onr1, onr1) end -""" - G2MomCorrelator(d::Int,c=:cross) <: AbstractOperator{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: -```math -\\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: -```math -\\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 -```math -\\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 - -* [`BoseHubbardMom1D2C`](@ref) -* [`BoseFS2C`](@ref) -* [`G2RealCorrelator`](@ref) -* [`G2RealSpace`](@ref) -* [`AbstractOperator`](@ref) -* [`AllOverlaps`](@ref) -""" -struct G2MomCorrelator{C} <: AbstractOperator{ComplexF64} - d::Int - - function G2MomCorrelator(d,c=:cross) - if c == :first - return new{1}(d) - elseif c == :second - return new{2}(d) - elseif c == :cross - return new{3}(d) - else - throw(ArgumentError("Unknown instruction for G2MomCorrelator!")) - end - end -end -@deprecate G2Correlator(d,c) G2MomCorrelator(d,c) -@deprecate G2Correlator(d) G2MomCorrelator(d) - -function Interfaces.allows_address_type(g2m::G2MomCorrelator, ::Type{A}) where {A} - return num_modes(A) > g2m.d -end - -function Base.show(io::IO, g::G2MomCorrelator{C}) where C - if C == 1 - print(io, "G2MomCorrelator($(g.d),:first)") - elseif C == 2 - print(io, "G2MomCorrelator($(g.d),:second)") - elseif C == 3 - print(io, "G2MomCorrelator($(g.d),:cross)") - end -end - - -num_offdiagonals(g::G2MomCorrelator{1},addr::BoseFS2C) = num_offdiagonals(g, addr.bsa) -num_offdiagonals(g::G2MomCorrelator{2},addr::BoseFS2C) = num_offdiagonals(g, addr.bsb) - -function num_offdiagonals(g::G2MomCorrelator{3}, addr::BoseFS2C) - m = num_modes(addr) - sa = num_occupied_modes(addr.bsa) - sb = num_occupied_modes(addr.bsb) - return sa*(m-1)*sb -end - -function num_offdiagonals(g::G2MomCorrelator, addr::SingleComponentFockAddress) - m = num_modes(addr) - singlies, doublies = num_singly_doubly_occupied_sites(addr) - return singlies*(singlies-1)*(m - 2) + doublies*(m - 1) -end - -diagonal_element(g::G2MomCorrelator{1}, addr::BoseFS2C) = diagonal_element(g, addr.bsa) -diagonal_element(g::G2MomCorrelator{2}, addr::BoseFS2C) = diagonal_element(g, addr.bsb) - -function diagonal_element(g::G2MomCorrelator{3}, addr::BoseFS2C{NA,NB,M,AA,AB}) where {NA,NB,M,AA,AB} - onrep_a = onr(addr.bsa) - onrep_b = onr(addr.bsb) - gd = 0 - for p in 1:M - iszero(onrep_b[p]) && continue - for k in 1:M - gd += onrep_a[k]*onrep_b[p] # b†_p b_p a†_k a_k - end - end - return ComplexF64(gd/M) -end - -function diagonal_element(g::G2MomCorrelator, addr::SingleComponentFockAddress) - M = num_modes(addr) - onrep = onr(addr) - gd = 0 - for p in 1:M - iszero(onrep[p]) && continue - for k in 1:M - gd += onrep[k]*onrep[p] # a†_p a_p a†_k a_k - end - end - return ComplexF64(gd/M) -end - -function get_offdiagonal(g::G2MomCorrelator{1}, addr::A, chosen)::Tuple{A,ComplexF64} where A<:BoseFS2C - new_bsa, elem = get_offdiagonal(g, addr.bsa, chosen) - return A(new_bsa,addr.bsb), elem -end - -function get_offdiagonal(g::G2MomCorrelator{2}, addr::A, chosen)::Tuple{A,ComplexF64} where A<:BoseFS2C - new_bsb, elem = get_offdiagonal(g, addr.bsb, chosen) - return A(addr.bsa, new_bsb), elem -end - -function get_offdiagonal( - g::G2MomCorrelator{3}, - addr::A, - chosen, - ma=OccupiedModeMap(addr.bsa), - mb=OccupiedModeMap(addr.bsb), -)::Tuple{A,ComplexF64} where {A<:BoseFS2C} - - m = num_modes(addr) - new_bsa, new_bsb, gamma, _, _, Δp = momentum_transfer_excitation( - addr.bsa, addr.bsb, chosen, ma, mb - ) - gd = exp(-im*g.d*Δp*2π/m)*gamma - return A(new_bsa, new_bsb), ComplexF64(gd/m) -end - -function get_offdiagonal( - g::G2MomCorrelator, - addr::A, - chosen, -)::Tuple{A,ComplexF64} where {A<:SingleComponentFockAddress} - M = num_modes(addr) - new_add, gamma, Δp = momentum_transfer_excitation(addr, chosen, OccupiedModeMap(addr)) - gd = exp(-im*g.d*Δp*2π/M)*gamma - return new_add, ComplexF64(gd/M) -end - """ SuperfluidCorrelator(d::Int) <: AbstractOperator{Float64} diff --git a/src/Hamiltonians/geometry.jl b/src/Hamiltonians/geometry.jl index 5d2992c4c..cd2f45db8 100644 --- a/src/Hamiltonians/geometry.jl +++ b/src/Hamiltonians/geometry.jl @@ -2,7 +2,8 @@ CubicGrid(dims::NTuple{D,Int}, fold::NTuple{D,Bool}) Represents a `D`-dimensional grid. Used to define a cubic lattice and boundary conditions -for some [`AbstractHamiltonian`](@ref)s. The type instance can be used to convert between +for some [`AbstractHamiltonian`](@ref)s, e.g. with the keyword argument `geometry` when +constructing a [`HubbardRealSpace`](@ref). The type instance can be used to convert between cartesian vector indices (tuples or `SVector`s) and linear indices (integers). When indexed with vectors, it folds them back into the grid if the out-of-bounds dimension is periodic and 0 otherwise (see example below). @@ -38,7 +39,8 @@ julia> geo[(3,4)] # returns 0 if out of bounds ``` See also [`PeriodicBoundaries`](@ref), [`HardwallBoundaries`](@ref) and -[`LadderBoundaries`](@ref) for special-case constructors. +[`LadderBoundaries`](@ref) for special-case constructors. See also +[`HubbardRealSpace`](@ref) and [`G2RealSpace`](@ref). """ struct CubicGrid{D,Dims,Fold} function CubicGrid( @@ -56,7 +58,7 @@ CubicGrid(args::Vararg{Int}) = CubicGrid(args) PeriodicBoundaries(dims...) -> CubicGrid PeriodicBoundaries(dims) -> CubicGrid -Return `CubicGrid` with all dimensions periodic. Equivalent to `CubicGrid(dims)`. +Return a [`CubicGrid`](@ref) with all dimensions periodic. Equivalent to `CubicGrid(dims)`. """ function PeriodicBoundaries(dims::NTuple{D,Int}) where {D} return CubicGrid(dims, ntuple(Returns(true), Val(D))) @@ -67,7 +69,7 @@ PeriodicBoundaries(dims::Vararg{Int}) = PeriodicBoundaries(dims) HardwallBoundaries(dims...) -> CubicGrid HardwallBoundaries(dims) -> CubicGrid -Return `CubicGrid` with all dimensions non-periodic. Equivalent to +Return a [`CubicGrid`](@ref) with all dimensions non-periodic. Equivalent to `CubicGrid(dims, (false, false, ...))`. """ function HardwallBoundaries(dims::NTuple{D,Int}) where {D} @@ -79,8 +81,8 @@ HardwallBoundaries(dims::Vararg{Int}) = HardwallBoundaries(dims) LadderBoundaries(dims...) -> CubicGrid LadderBoundaries(dims) -> CubicGrid -Return `CubicGrid` where the first dimension is dimensions non-periodic and the rest are -periodic. Equivalent to `CubicGrid(dims, (true, false, ...))`. +Return a [`CubicGrid`](@ref) where the first dimension is dimensions non-periodic and the +rest are periodic. Equivalent to `CubicGrid(dims, (true, false, ...))`. """ function LadderBoundaries(dims::NTuple{D,Int}) where {D} return CubicGrid(dims, ntuple(>(1), Val(D))) @@ -100,13 +102,15 @@ fold(g::CubicGrid{<:Any,<:Any,Fold}) where {Fold} = Fold num_dimensions(geom::LatticeCubicGrid) Return the number of dimensions of the lattice in this geometry. + +See also [`CubicGrid`](@ref). """ num_dimensions(::CubicGrid{D}) where {D} = D """ fold_vec(g::CubicGrid{D}, vec::SVector{D,Int}) -> SVector{D,Int} -Use the CubicGrid to fold the `vec` in each dimension. If folding is disabled in a +Use the [`CubicGrid`](@ref) to fold the `vec` in each dimension. If folding is disabled in a dimension, and the vector is allowed to go out of bounds. ```julia @@ -179,8 +183,8 @@ end """ Displacements(geometry::CubicGrid) <: AbstractVector{SVector{D,Int}} -Return all valid offset vectors in a [`CubicGrid`](@ref). If `center=true` the (0,0) displacement is -placed at the centre of the array. +Return all valid offset vectors in a [`CubicGrid`](@ref). If `center=true` the (0,0) +displacement is placed at the centre of the array. ```jldoctest; setup=:(using Rimu.Hamiltonians: Displacements) julia> geometry = CubicGrid((3,4)); @@ -222,6 +226,8 @@ end neighbor_site(geom::CubicGrid, site, i) Find the `i`-th neighbor of `site` in the geometry. If the move is illegal, return 0. + +See also [`CubicGrid`](@ref). """ function neighbor_site(g::CubicGrid{D}, mode, chosen) where {D} # TODO: reintroduce LadderBoundaries small dimensions diff --git a/src/Hamiltonians/ho-cart-tools.jl b/src/Hamiltonians/ho-cart-tools.jl index 9ca720829..ecb58529f 100644 --- a/src/Hamiltonians/ho-cart-tools.jl +++ b/src/Hamiltonians/ho-cart-tools.jl @@ -1,18 +1,18 @@ """ - get_all_blocks(h::Union{HOCartesianContactInteractions,HOCartesianEnergyConservedPerDim}; - target_energy = nothing, - max_energy = nothing, - max_blocks = nothing, + 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 +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`](@ref) -* `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 +* `addr`: first address that generates the block with e.g. [`BasisSetRepresentation`](@ref) +* `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 `NamedTuple`s. * `t_basis`: time to generate the basis for each block @@ -21,22 +21,22 @@ 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 +* `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 +* `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. +* 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`](@ref) with an appropriate `filter` to generate the blocks. +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 +[`BasisSetRepresentation`](@ref) with an appropriate `filter` to generate the blocks. """ -function get_all_blocks(h::Union{HOCartesianContactInteractions{D,<:Any,B},HOCartesianEnergyConservedPerDim{D}}; - target_energy = nothing, - max_energy = nothing, - max_blocks = nothing, +function get_all_blocks(h::Union{HOCartesianContactInteractions{D,<:Any,B},HOCartesianEnergyConservedPerDim{D}}; + target_energy = nothing, + max_energy = nothing, + max_blocks = nothing, method = :vertices, kwargs... ) where {D,B} @@ -63,7 +63,7 @@ function get_all_blocks(h::Union{HOCartesianContactInteractions{D,<:Any,B},HOCar (!isnothing(max_energy) && max_energy > E0 + M) || (!isnothing(target_energy) && target_energy > E0 + M) throw(ArgumentError("requested energy range not accessible by given Hamiltonian")) - end + end end if method == :vertices @@ -75,10 +75,10 @@ function get_all_blocks(h::Union{HOCartesianContactInteractions{D,<:Any,B},HOCar end end -function get_all_blocks_vertices(h; - target_energy = nothing, - max_energy = nothing, - max_blocks = nothing, +function get_all_blocks_vertices(h; + target_energy = nothing, + max_energy = nothing, + max_blocks = nothing, save_to_file = nothing, kwargs... ) @@ -124,10 +124,10 @@ function get_all_blocks_vertices(h; end # old version - issues with GC due to allocating many small vectors -function get_all_blocks_comb(h; - target_energy = nothing, - max_energy = nothing, - max_blocks = nothing, +function get_all_blocks_comb(h; + target_energy = nothing, + max_energy = nothing, + max_blocks = nothing, save_to_file = nothing, kwargs... ) @@ -173,9 +173,9 @@ end """ parity_block_seed_addresses(h::HOCartesianContactInteractions{D}) -Get a vector of addresses that each have different parity with respect to +Get a vector of addresses that each have different parity with respect to the trap geometry defined by the Hamiltonian `H`. The result will have `2^D` -`BoseFS` addresses for a `D`-dimensional trap. This is useful for +`BoseFS` addresses for a `D`-dimensional trap. This is useful for [`HOCartesianContactInteractions`](@ref) with option `block_by_level = false`. """ function parity_block_seed_addresses(h::HOCartesianContactInteractions{D,A}) where {D,A} @@ -189,7 +189,7 @@ function parity_block_seed_addresses(h::HOCartesianContactInteractions{D,A}) whe BoseFS(P, index => 1, 1 => N-1) ) end - + return seeds end @@ -197,7 +197,8 @@ end function get_all_blocks_parity(h::HOCartesianContactInteractions{D,<:Any,B}) where {D,B} # check if H blocks by level B && throw(ArgumentError("use `get_all_blocks` instead")) - bs_estimate = prod(h.S) # rough upper bound assuming `BasisSetRep` will filter by energy cutoff + bs_estimate = prod(h.S) + # rough upper bound assuming `BasisSetRepresentation` will filter by energy cutoff df = DataFrame() for (block_id, addr) in enumerate(parity_block_seed_addresses(h)) t = vcat(map(p -> [p.mode for _ in 1:p.occnum], OccupiedModeMap(addr))...) @@ -210,23 +211,23 @@ end """ 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 +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`. """ function fock_to_cart( - addr::SingleComponentFockAddress{N}, - S::NTuple{D,Int}; + addr::SingleComponentFockAddress{N}, + S::NTuple{D,Int}; zero_index = true ) where {N,D} prod(S) == num_modes(addr) || throw(ArgumentError("Specified cartesian states are incompatible with address")) states = CartesianIndices(S) cart = vcat(map( - p -> [Tuple(states[p.mode]) .- Int(zero_index) for _ in 1:p.occnum], + p -> [Tuple(states[p.mode]) .- Int(zero_index) for _ in 1:p.occnum], OccupiedModeMap(addr))...) return SVector{N,NTuple{D,Int}}(cart) -end \ No newline at end of file +end diff --git a/src/InterfaceTests.jl b/src/InterfaceTests.jl new file mode 100644 index 000000000..44ca7c061 --- /dev/null +++ b/src/InterfaceTests.jl @@ -0,0 +1,285 @@ +""" +The module `Rimu.InterfaceTests` provides functions to test compliance with the +[`AbstractObservable`](@ref), [`AbstractOperator`](@ref), and [`AbstractHamiltonian`](@ref) +interfaces. Load the module with `using Rimu.InterfaceTests`. + +The module exports the following functions: +- [`test_observable_interface`](@ref Rimu.InterfaceTests.test_observable_interface) +- [`test_operator_interface`](@ref Rimu.InterfaceTests.test_operator_interface) +- [`test_hamiltonian_interface`](@ref Rimu.InterfaceTests.test_hamiltonian_interface) +- [`test_hamiltonian_structure`](@ref Rimu.InterfaceTests.test_hamiltonian_structure) +""" +module InterfaceTests + +using Test: Test, @test, @testset, @test_throws +using Rimu: Rimu, DVec, Interfaces, LOStructure, IsHermitian, IsDiagonal, AdjointKnown, + Hamiltonians, num_offdiagonals, allows_address_type, offdiagonals, random_offdiagonal, + diagonal_element, dimension, dot_from_right, IsDeterministic, starting_address, PDVec, + sparse, scale!, scalartype +using Rimu.Hamiltonians: AbstractHamiltonian, AbstractOperator, AbstractObservable, + AbstractOffdiagonals +using LinearAlgebra: dot, mul!, isdiag, ishermitian + +export test_observable_interface, test_operator_interface, test_hamiltonian_interface, + test_hamiltonian_structure + +""" + test_observable_interface(obs, addr) + +This function tests compliance with the [`AbstractObservable`](@ref) interface for an +observable `obs` at address `addr` (typically +[`<: AbstractFockAddress`](@ref Rimu.BitStringAddresses.AbstractFockAddress)) by checking +that all required methods are defined. + +The following properties are tested: +- `dot(v, obs, v)` returns a value of the same type as the `eltype` of the observable +- `LOStructure` is set consistently + +### Example +```julia-doctest +julia> using Rimu.InterfaceTests + +julia> test_observable_interface(ReducedDensityMatrix(2), FermiFS(1,0,1,1)); +Test Summary: | Pass Total Time +Observable interface: ReducedDensityMatrix | 4 4 0.0s +``` + +See also [`AbstractObservable`](@ref), [`test_operator_interface`](@ref), +[`test_hamiltonian_interface`](@ref). +""" +function test_observable_interface(obs, addr) + @testset "`AbstractObservable` interface test: $(nameof(typeof(obs)))" begin + @testset "three way dot" begin # this works with vector valued operators + v = DVec(addr => scalartype(obs)(2)) + @test dot(v, obs, v) isa eltype(obs) + @test dot(v, obs, v) ≈ Interfaces.dot_from_right(v, obs, v) + end + @testset "LOStructure" begin + @test LOStructure(obs) isa LOStructure + if LOStructure(obs) isa IsHermitian + @test obs' === obs + elseif LOStructure(obs) isa IsDiagonal + @test num_offdiagonals(obs, addr) == 0 + if scalartype(obs) <: Real + @test obs' === obs + end + elseif LOStructure(obs) isa AdjointKnown + @test begin + obs' + true + end # make sure no error is thrown + else + @test_throws ArgumentError obs' + end + end + end +end + +""" + test_operator_interface(op, addr; test_spawning=true) + +This function tests compliance with the [`AbstractOperator`](@ref) interface for an operator +`op` at address `addr` (typically +[`<: AbstractFockAddress`](@ref Rimu.BitStringAddresses.AbstractFockAddress)) by +checking that all required methods are defined. + +If `test_spawning` is `true`, tests are performed that require `offdiagonals` to return an +`Hamiltonians.AbstractOffDiagonals`, which is a prerequisite for using the `spawn!` +function. Otherwise, the spawning tests are skipped. + +The following properties are tested: +- `diagonal_element` returns a value of the same type as the `eltype` of the operator +- `offdiagonals` behaves like an `AbstractVector` +- `num_offdiagonals` returns the correct number of offdiagonals +- `random_offdiagonal` returns a tuple with the correct types +- `mul!` and `dot` work as expected +- `dimension` returns a consistent value +- the [`AbstractObservable`](@ref) interface is tested + +### Example +```julia-doctest +julia> using Rimu.InterfaceTests + +julia> test_operator_interface(SuperfluidCorrelator(3), BoseFS(1, 2, 3, 1)); +Test Summary: | Pass Total Time +Observable interface: SuperfluidCorrelator | 4 4 0.0s +Test Summary: | Pass Total Time +allows_address_type | 1 1 0.0s +Test Summary: | Pass Total Time +Operator interface: SuperfluidCorrelator | 9 9 0.0s +``` + +See also [`AbstractOperator`](@ref), [`test_observable_interface`](@ref), +[`test_hamiltonian_interface`](@ref). +""" +function test_operator_interface(op, addr; test_spawning=true) + test_observable_interface(op, addr) + @testset "`AbstractOperator` interface test: $(nameof(typeof(op)))" begin + + @testset "allows_address_type" begin + @test allows_address_type(op, addr) + end + @testset "Operator interface: $(nameof(typeof(op)))" begin + @testset "diagonal_element" begin + @test diagonal_element(op, addr) isa eltype(op) + @test eltype(diagonal_element(op, addr)) == scalartype(op) + end + @testset "offdiagonals" begin + # `get_offdiagonal` is not mandatory and thus not tested + ods = offdiagonals(op, addr) + vec_ods = collect(ods) + eltype(vec_ods) == Tuple{typeof(addr),eltype(op)} == eltype(ods) + @test length(vec_ods) ≤ num_offdiagonals(op, addr) + end + if test_spawning + @testset "spawning" begin + ods = offdiagonals(op, addr) + @test ods isa AbstractOffdiagonals{typeof(addr),eltype(op)} + @test ods isa AbstractVector + @test size(ods) == (num_offdiagonals(op, addr),) + if length(ods) > 0 + @test random_offdiagonal(op, addr) isa Tuple{typeof(addr),<:Real,eltype(op)} + end + end + end + @testset "mul!" begin # this works with vector valued operators + v = DVec(addr => scalartype(op)(2)) + w = empty(v, eltype(op); style=IsDeterministic{scalartype(op)}()) + mul!(w, op, v) # operator vector product + @test dot(v, op, v) ≈ Interfaces.dot_from_right(v, op, v) ≈ dot(v, w) + end + @testset "dimension" begin + @test dimension(addr) ≥ dimension(op, addr) + end + end + end +end + +""" + test_hamiltonian_interface(h, addr=starting_address(h); test_spawning=true) + +The main purpose of this test function is to check that all required methods of the +[`AbstractHamiltonian`](@ref) interface are defined and work as expected. + +Set `test_spawning=false` to skip tests that require [`offdiagonals`](@ref) to return an +`AbstractVector`. + +This function also tests the following properties of the Hamiltonian: +- `dimension(h) ≥ dimension(h, addr)` +- `scalartype(h) === eltype(h)` +- Hamiltonian action on a vector <: `AbstractDVec` +- `starting_address` returns an [`allows_address_type`](@ref) address +- `LOStructure` is one of `IsDiagonal`, `IsHermitian`, `AdjointKnown` +- the [`AbstractOperator`](@ref) interface is tested +- the [`AbstractObservable`](@ref) interface is tested + +### Example +```julia-doctest +julia> using Rimu.InterfaceTests + +julia> test_hamiltonian_interface(HubbardRealSpace(BoseFS(2,0,3,1))); +Test Summary: | Pass Total Time +Observable interface: HubbardRealSpace | 4 4 0.0s +Test Summary: | Pass Total Time +allows_address_type | 1 1 0.0s +Test Summary: | Pass Total Time +Operator interface: HubbardRealSpace | 9 9 0.0s +Test Summary: | Pass Total Time +allows_address_type | 1 1 0.0s +Test Summary: | Pass Total Time +Hamiltonians-only tests with HubbardRealSpace | 6 6 0.0s +``` + +See also [`test_operator_interface`](@ref), [`test_observable_interface`](@ref). +""" +function test_hamiltonian_interface(h, addr=starting_address(h); test_spawning=true) + test_operator_interface(h, addr; test_spawning) + @testset "`AbstractHamiltonian` interface test: $(nameof(typeof(h)))" begin + + @testset "allows_address_type on starting_address" begin + @test allows_address_type(h, starting_address(h)) + end + @testset "Hamiltonians-only tests with $(nameof(typeof(h)))" begin + # starting_address is specific to Hamiltonians + @test allows_address_type(h, starting_address(h)) + + @test dimension(h) ≥ dimension(h, addr) + + # Hamiltonians can only have scalar eltype + @test scalartype(h) === eltype(h) + + # Hamiltonian action on a vector + v = DVec(addr => scalartype(h)(2)) + v1 = similar(v) + mul!(v1, h, v) + v2 = h * v + v3 = similar(v) + h(v3, v) + v4 = h(v) + @test v1 == v2 == v3 == v4 + v5 = DVec(addr => diagonal_element(h, addr)) + for (addr, val) in offdiagonals(h, addr) + v5[addr] += val + end + scale!(v5, scalartype(h)(2)) + v5[addr] = v5[addr] # remove possible 0.0 from the diagonal + @test v5 == v1 + + if test_spawning && scalartype(h) <: Real + # applying an operator on a PDVec uses spawn!, which requires + # offdiagonals to be an AbstractVector + # currently this only works for real operators as spawn! is not + # implemented for complex operators + pv = PDVec(addr => scalartype(h)(2)) + pv1 = h(pv) + @test dot(pv1, h, pv) ≈ Interfaces.dot_from_right(pv1, h, pv) ≈ dot(v1, v1) + end + end + end +end + +""" + test_hamiltonian_structure(h::AbstractHamiltonian; sizelim=20) + +Test the `LOStructure` of a small Hamiltonian `h` by instantiating it as a sparse matrix and +checking whether the structure of the matrix is constistent with the result of +`LOStructure(h)` and the `eltype` is consistent with `eltype(h)`. + +This function is intended to be used in automated test for small Hamiltonians where +instantiating the matrix is quick. A warning will print if the dimension of the Hamiltonian +is larger than `20`. + +### Example +```julia-doctest +julia> using Rimu.InterfaceTests + +julia> test_hamiltonian_structure(HubbardRealSpace(BoseFS(2,0,1))); +Test Summary: | Pass Total Time +structure | 4 4 0.0s +``` +""" +function test_hamiltonian_structure(h::AbstractHamiltonian; sizelim=20) + @testset "Hamiltonian structure: $(nameof(typeof(h)))" begin + d = dimension(h) + d > 20 && @warn "This function is intended for small Hamiltonians. The dimension is $d." + m = sparse(h; sizelim) + @test eltype(m) === eltype(h) + if LOStructure(h) == IsDiagonal() + @test isdiag(m) + elseif LOStructure(h) == IsHermitian() + @test h' == h + @test h' === h + @test ishermitian(m) + elseif LOStructure(h) == AdjointKnown() + @test m' == sparse(h') + end + if !ishermitian(m) + @test LOStructure(h) != IsHermitian() + if LOStructure(h) == IsDiagonal() + @test !isreal(m) + @test !(eltype(h) <: Real) + end + end + end +end +end diff --git a/src/Rimu.jl b/src/Rimu.jl index ae40a363b..7a190185c 100644 --- a/src/Rimu.jl +++ b/src/Rimu.jl @@ -104,4 +104,6 @@ include("pmc_simulation.jl") include("lomc.jl") # top level +include("InterfaceTests.jl") + end # module diff --git a/src/strategies_and_params/poststepstrategy.jl b/src/strategies_and_params/poststepstrategy.jl index 49dfb282b..1bdd33240 100644 --- a/src/strategies_and_params/poststepstrategy.jl +++ b/src/strategies_and_params/poststepstrategy.jl @@ -286,7 +286,7 @@ end function single_particle_density(add::SingleComponentFockAddress; component=0) return float.(Tuple(onr(add))) end -function single_particle_density(add::Union{CompositeFS,BoseFS2C}; component=0) +function single_particle_density(add::Union{CompositeFS}; component=0) if component == 0 return float.(Tuple(sum(onr(add)))) else diff --git a/test/BitStringAddresses.jl b/test/BitStringAddresses.jl index 1d7480e01..efe35896d 100644 --- a/test/BitStringAddresses.jl +++ b/test/BitStringAddresses.jl @@ -450,24 +450,6 @@ end @test num_particles(FermiFS2C(3, 1 => 1)) == 1 @test num_particles(FermiFS2C(3, 1 => -1)) == 1 end - - @testset "BoseFS2C" begin - fs1 = BoseFS2C((1,1,1,0,0,0,0), (1,1,1,0,5,0,0)) - @test num_modes(fs1) == 7 - @test num_components(fs1) == 2 - @test num_particles(fs1) == 11 - @test eval(Meta.parse(repr(fs1))) == fs1 - @test BoseFS2C(parse_address(sprint(show, fs1; context=:compact => true))) == fs1 - @test onr(fs1) == ([1, 1, 1, 0, 0, 0, 0], [1, 1, 1, 0, 5, 0, 0]) - - fs2 = BoseFS2C(BoseFS((0,0,0,0,0,0,3)), BoseFS((0,2,1,0,5,0,0))) - @test fs1 < fs2 - - @test_throws MethodError BoseFS2C(BoseFS((1,1)), BoseFS((1,1,1))) - @test BoseFS2C(CompositeFS(BoseFS((1,2)), BoseFS((3,1)))) == BoseFS2C((1,2), (3,1)) - @test CompositeFS(BoseFS2C(BoseFS((1,2)), BoseFS((3,1)))) == - CompositeFS(BoseFS((1,2)), BoseFS((3,1))) - end end @testset "OccupationNumberFS functions" begin diff --git a/test/Hamiltonians.jl b/test/Hamiltonians.jl index caee1a335..03451a2b2 100644 --- a/test/Hamiltonians.jl +++ b/test/Hamiltonians.jl @@ -7,6 +7,8 @@ using DataFrames using Suppressor using StaticArrays using Rimu.Hamiltonians: TransformUndoer, AbstractOffdiagonals +using Rimu.InterfaceTests: test_observable_interface, test_operator_interface, + test_hamiltonian_interface, test_hamiltonian_structure function exact_energy(ham) dv = DVec(starting_address(ham) => 1.0) @@ -14,175 +16,6 @@ function exact_energy(ham) return all_results[1][1] end -""" - test_observable_interface(obs, addr) - -This function tests the interface of an observable `obs` at address `addr` by checking that -all required methods are defined. -""" -function test_observable_interface(obs, addr) - @testset "Observable interface: $(nameof(typeof(obs)))" begin - @testset "three way dot" begin # this works with vector valued operators - v = DVec(addr => scalartype(obs)(2)) - @test dot(v, obs, v) isa eltype(obs) - @test dot(v, obs, v) ≈ Interfaces.dot_from_right(v, obs, v) - end - @testset "LOStructure" begin - @test LOStructure(obs) isa LOStructure - if LOStructure(obs) isa IsHermitian - @test obs' === obs - elseif LOStructure(obs) isa IsDiagonal - @test num_offdiagonals(obs, addr) == 0 - if scalartype(obs) <: Real - @test obs' === obs - end - elseif LOStructure(obs) isa AdjointKnown - @test begin obs'; true; end # make sure no error is thrown - else - @test_throws ArgumentError obs' - end - end - @testset "show" begin - # Check that the result of show can be pasted into the REPL - @test eval(Meta.parse(repr(obs))) == obs - end - end -end - -""" - test_operator_interface(op, addr; test_spawning=true) - -This function tests the interface of an operator `op` at address `addr` by checking that all -required methods are defined. - -If `test_spawning` is `true`, tests are performed that require `offdiagonals` to return an -`Hamiltonians.AbstractOffDiagonals`, which is a prerequisite for using the `spawn!` -function. Otherwise, the spawning tests are skipped. -""" -function test_operator_interface(op, addr; test_spawning=true) - test_observable_interface(op, addr) - - @testset "allows_address_type" begin - @test allows_address_type(op, addr) - end - @testset "Operator interface: $(nameof(typeof(op)))" begin - @testset "diagonal_element" begin - @test diagonal_element(op, addr) isa eltype(op) - @test eltype(diagonal_element(op, addr)) == scalartype(op) - end - @testset "offdiagonals" begin - # `get_offdiagonal` is not mandatory and thus not tested - ods = offdiagonals(op, addr) - vec_ods = collect(ods) - eltype(vec_ods) == Tuple{typeof(addr), eltype(op)} == eltype(ods) - @test length(vec_ods) ≤ num_offdiagonals(op, addr) - end - if test_spawning - @testset "spawning" begin - ods = offdiagonals(op, addr) - @test ods isa AbstractOffdiagonals{typeof(addr),eltype(op)} - @test ods isa AbstractVector - @test size(ods) == (num_offdiagonals(op, addr),) - if length(ods) > 0 - @test random_offdiagonal(op, addr) isa Tuple{typeof(addr), <:Real, eltype(op)} - end - end - end - @testset "mul!" begin # this works with vector valued operators - v = DVec(addr => scalartype(op)(2)) - w = empty(v, eltype(op); style=IsDeterministic{scalartype(op)}()) - mul!(w, op, v) # operator vector product - @test dot(v, op, v) ≈ Interfaces.dot_from_right(v, op, v) ≈ dot(v, w) - end - @testset "dimension" begin - @test dimension(addr) ≥ dimension(op, addr) - end - end -end - -""" - test_hamiltonian_interface(H, addr=starting_address(H)) - -The main purpose of this test function is to check that all required methods are defined. -""" -function test_hamiltonian_interface(H, addr=starting_address(H); test_spawning=true) - test_operator_interface(H, addr; test_spawning) - - @testset "allows_address_type" begin - @test allows_address_type(H, addr) - end - @testset "Hamiltonians-only tests with $(nameof(typeof(H)))" begin - # starting_address is specific to Hamiltonians - @test allows_address_type(H, starting_address(H)) - - @test dimension(H) ≥ dimension(H, addr) - - # Hamiltonians can only have scalar eltype - @test scalartype(H) === eltype(H) - - # Hamiltonian action on a vector - v = DVec(addr => scalartype(H)(2)) - v1 = similar(v) - mul!(v1, H, v) - v2 = H * v - v3 = similar(v) - H(v3, v) - v4 = H(v) - @test v1 == v2 == v3 == v4 - v5 = DVec(addr => diagonal_element(H, addr)) - for (addr, val) in offdiagonals(H, addr) - v5[addr] += val - end - scale!(v5, scalartype(H)(2)) - v5[addr] = v5[addr] # remove possible 0.0 from the diagonal - @test v5 == v1 - - if test_spawning && scalartype(H) <: Real - # applying an operator on a PDVec uses spawn!, which requires - # offdiagonals to be an AbstractVector - # currently this only works for real operators as spawn! is not - # implemented for complex operators - pv = PDVec(addr => scalartype(H)(2)) - pv1 = H(pv) - @test dot(pv1, H, pv) ≈ Interfaces.dot_from_right(pv1, H, pv) ≈ dot(v1, v1) - end - - end -end - -""" - test_hamiltonian_structure(H) - -Test the `LOStructure` of a small Hamiltonian `H` by instantiating it as a sparse matrix and -checking whether the structure of the matrix is constistent with the result of -`LOStructure(H)` and the `eltype` is consistent with `eltype(H)`. - -This function is intended to be used for small Hamiltonians where instantiating the matrix -is quick. -""" -function test_hamiltonian_structure(H) - d = dimension(H) - d > 20 && @warn "This function is intended for small Hamiltonians. The dimension is $d." - m = sparse(H) - @test eltype(m) === eltype(H) - if LOStructure(H) == IsDiagonal() - @test isdiag(m) - elseif LOStructure(H) == IsHermitian() - @test H' == H - @test H' === H - @test ishermitian(m) - elseif LOStructure(H) == AdjointKnown() - @test m' == sparse(H') - end - if !ishermitian(m) - @test LOStructure(H) != IsHermitian() - if LOStructure(H) == IsDiagonal() - @test !isreal(m) - @test !(eltype(H) <: Real) - end - end - nothing -end @testset "Hamiltonian interface tests" begin for H in ( HubbardReal1D(BoseFS((1, 2, 3, 4)); u=1.0, t=2.0), @@ -205,10 +38,8 @@ end FermiFS((1, 1, 1, 1, 1, 0, 0, 0)), FermiFS((1, 1, 1, 1, 0, 0, 0, 0)), ); t=[1, 2], u=[0 3; 3 0] - ), BoseHubbardReal1D2C(BoseFS2C((1, 2, 3), (1, 0, 0))), - BoseHubbardMom1D2C(BoseFS2C((1, 2, 3), (1, 0, 0))), + ), GutzwillerSampling(HubbardReal1D(BoseFS((1, 2, 3)); u=6); g=0.3), - GutzwillerSampling(BoseHubbardMom1D2C(BoseFS2C((3, 2, 1), (1, 2, 3)); ua=6); g=0.3), GutzwillerSampling(HubbardReal1D(BoseFS((1, 2, 3)); u=6 + 2im); g=0.3), MatrixHamiltonian(Float64[1 2; 2 0]), GutzwillerSampling(MatrixHamiltonian([1.0 2.0; 2.0 0.0]); g=0.3), @@ -221,7 +52,6 @@ end HubbardMom1DEP(CompositeFS(FermiFS((0, 1, 1, 0, 0)), FermiFS((0, 0, 1, 0, 0))), v_ho=5), ParitySymmetry(HubbardRealSpace(CompositeFS(BoseFS((1, 2, 0)), FermiFS((0, 1, 0))))), TimeReversalSymmetry(HubbardMom1D(FermiFS2C((1, 0, 1), (0, 1, 1)))), - TimeReversalSymmetry(BoseHubbardMom1D2C(BoseFS2C((0, 1, 1), (1, 0, 1)))), Stoquastic(HubbardMom1D(BoseFS((0, 5, 0)))), momentum(HubbardMom1D(BoseFS((0, 5, 0)))), HOCartesianContactInteractions(BoseFS((2, 0, 0, 0))), @@ -233,6 +63,9 @@ end ) # test_hamiltonian_interface(H; test_spawning=false) test_hamiltonian_interface(H; test_spawning=!(H isa HOCartesianContactInteractions)) + # Check that the result of show can be pasted into the REPL + @test eval(Meta.parse(repr(H))) == H + end end @@ -246,15 +79,17 @@ end (ParticleNumberOperator(), BoseFS(1, 2, 3)), (ParticleNumberOperator(), FermiFS2C((1, 0, 1), (0, 1, 1))), (G2RealCorrelator(3), FermiFS2C((1, 0, 1, 1), (0, 1, 1, 0))), - (G2MomCorrelator(3), BoseFS(1, 2, 0, 3, 0, 4, 0, 1)), (SuperfluidCorrelator(3), BoseFS(1, 2, 3, 1)), (StringCorrelator(3), BoseFS(1, 0, 3, 1)), (DensityMatrixDiagonal(3), FermiFS(1, 0, 1)), (SingleParticleExcitation(2, 3), BoseFS(1, 2, 3, 4)), (TwoParticleExcitation(3, 2, 1, 4), BoseFS(1, 2, 3, 4)), (Momentum(1), BoseFS(1, 2, 3, 4)), + (G2MomCorrelator(3), BoseFS(1, 2, 0, 3, 0, 4, 0, 1)) ] test_operator_interface(op, addr) + # Check that the result of show can be pasted into the REPL + @test eval(Meta.parse(repr(op))) == op end end @@ -425,36 +260,6 @@ end end end -@testset "2C model properties" begin - flip(b) = BoseFS2C(b.bsb, b.bsa) - addr1 = near_uniform(BoseFS2C{1,100,20}) - addr2 = near_uniform(BoseFS2C{100,1,20}) - - for Hamiltonian in (BoseHubbardReal1D2C, BoseHubbardMom1D2C) - @testset "$Hamiltonian" begin - H1 = BoseHubbardReal1D2C(addr1; ta=1.0, tb=2.0, ua=0.5, ub=0.7, v=0.2) - H2 = BoseHubbardReal1D2C(addr2; ta=2.0, tb=1.0, ua=0.7, ub=0.5, v=0.2) - @test starting_address(H1) == addr1 - @test LOStructure(H1) == IsHermitian() - - hops1 = collect(offdiagonals(H1, addr1)) - hops2 = collect(offdiagonals(H2, addr2)) - sort!(hops1, by=a -> first(a).bsa) - sort!(hops2, by=a -> first(a).bsb) - - addrs1 = first.(hops1) - addrs2 = flip.(first.(hops2)) - values1 = last.(hops1) - values2 = last.(hops1) - @test addrs1 == addrs2 - @test values1 == values2 - - @test eval(Meta.parse(repr(H1))) == H1 - @test eval(Meta.parse(repr(H2))) == H2 - end - end -end - @testset "HubbardRealSpace" begin @testset "Constructor" begin bose = BoseFS((1, 2, 3, 4, 5, 6)) @@ -479,8 +284,6 @@ end comp; geometry=PeriodicBoundaries(3,2), v=[1 1; 1 1; 1 1], ) - @test_throws ArgumentError HubbardRealSpace(BoseFS2C((1,2,3), (3,2,1))) - @test_logs (:warn,) HubbardRealSpace(FermiFS((1,0)), u=[2]) @test_logs (:warn,) HubbardRealSpace( CompositeFS(BoseFS((1,1)), FermiFS((1,0))); u=[2 2; 2 2] @@ -523,12 +326,6 @@ end @test exact_energy(H1) == exact_energy(H2) end @testset "1D Bosons (2-component)" begin - add1 = BoseFS2C( - (1, 1, 1, 0, 0, 0), - (1, 0, 0, 0, 0, 0), - ) - H1 = BoseHubbardReal1D2C(add1, ua=2, v=3, tb=4) - add2 = CompositeFS( BoseFS((1, 1, 1, 0, 0, 0)), BoseFS((1, 0, 0, 0, 0, 0)), @@ -553,13 +350,11 @@ end ) H5 = HubbardRealSpace(add5, t=[4,1], u=[0 3; 3 2]) - E1 = exact_energy(H1) E2 = exact_energy(H2) E3 = exact_energy(H3) E4 = exact_energy(H4) E5 = exact_energy(H5) - @test E1 ≈ E2 rtol=0.0001 @test E2 ≈ E3 rtol=0.0001 @test E3 ≈ E4 rtol=0.0001 @test E4 ≈ E5 rtol=0.0001 @@ -698,7 +493,6 @@ end HubbardMom1D(BoseFS((2,2,2)), u=6), ExtendedHubbardReal1D(BoseFS((1,1,1,1,1,1,1,1,1,1,1,1)), u=6, t=2.0), ExtendedHubbardMom1D(BoseFS((1,1,1,1,1,1,1,1,1,1,1,1)), u=6, t=2.0), - BoseHubbardMom1D2C(BoseFS2C((1,2,3), (1,0,0)), ub=2.0), ) # GutzwillerSampling with parameter zero is exactly equal to the original H G = GutzwillerSampling(H, 0.0) @@ -729,7 +523,6 @@ end HubbardMom1D(BoseFS((2,2,2)), u=6), ExtendedHubbardReal1D(BoseFS((1,1,1,1,1,1,1,1,1,1,1,1)), u=6, t=2.0), ExtendedHubbardMom1D(BoseFS((1,1,1,1,1,1,1,1,1,1,1,1)), u=6, t=2.0), - # BoseHubbardMom1D2C(BoseFS2C((1,2,3), (1,0,0)), ub=2.0), # multicomponent not implemented for G2RealCorrelator ) # energy g = rand() @@ -753,10 +546,6 @@ end 0:m-1 ) @test all(g2vals ≈ g2transformed) - - # type promotion - G2mom = G2MomCorrelator(1) - @test eltype(Rimu.Hamiltonians.TransformUndoer(G, G2mom)) == eltype(G2mom) end end end @@ -813,7 +602,6 @@ end HubbardMom1D(BoseFS((2,2,2)), u=6), ExtendedHubbardReal1D(BoseFS((1,1,1,1,1,1,1,1,1,1,1,1)), u=6, t=2.0), ExtendedHubbardMom1D(BoseFS((1,1,1,1,1,1,1,1,1,1,1,1)), u=6, t=2.0), - # BoseHubbardMom1D2C(BoseFS2C((1,2,3), (1,0,0)), ub=2.0), # multicomponent not implemented for G2RealCorrelator ) # energy x = rand() @@ -833,10 +621,6 @@ end g2vals = map(d -> dot(dv, G2RealCorrelator(d), dv)/dot(dv, dv), 0:m-1) g2transformed = map(d -> dot(dv, Rimu.Hamiltonians.TransformUndoer(G,G2RealCorrelator(d)), dv)/dot(dv, fsq, dv), 0:m-1) @test all(g2vals ≈ g2transformed) - - # type promotion - G2mom = G2MomCorrelator(1) - @test eltype(Rimu.Hamiltonians.TransformUndoer(G, G2mom)) == eltype(G2mom) end end end @@ -883,7 +667,6 @@ end HubbardMom1D(BoseFS((2,2,2)), u=6), ExtendedHubbardReal1D(BoseFS((1,1,1,1,1,1,1,1,1,1,1,1)), u=6, t=2.0), ExtendedHubbardMom1D(BoseFS((1,1,1,1,1,1,1,1,1,1,1,1)), u=6, t=2.0), - BoseHubbardMom1D2C(BoseFS2C((1,2,3), (1,0,0)), ub=2.0), ) @test_throws ArgumentError Rimu.Hamiltonians.TransformUndoer(H) @test_throws ArgumentError Rimu.Hamiltonians.TransformUndoer(H, H) @@ -1096,55 +879,6 @@ using Rimu.Hamiltonians: circshift_dot end end - @testset "G2MomCorrelator" begin - # v0 is the exact ground state from BoseHubbardMom1D2C(aIni;ua=0,ub=0,v=0.1) - bfs1 = BoseFS([0, 2, 0]) - bfs2 = BoseFS([0, 1, 0]) - aIni = BoseFS2C(bfs1,bfs2) - v0 = DVec( - BoseFS2C((0, 2, 0), (0, 1, 0)) => 0.9999389545691221, - BoseFS2C((1, 1, 0), (0, 0, 1)) => -0.007812695959057453, - BoseFS2C((0, 1, 1), (1, 0, 0)) => -0.007812695959057453, - BoseFS2C((2, 0, 0), (1, 0, 0)) => 4.046694762039993e-5, - BoseFS2C((0, 0, 2), (0, 0, 1)) => 4.046694762039993e-5, - BoseFS2C((1, 0, 1), (0, 1, 0)) => 8.616127793651117e-5, - ) - g0 = G2MomCorrelator(0) - g1 = G2MomCorrelator(1) - g2 = G2MomCorrelator(2) - g3 = G2MomCorrelator(3) - @test imag(dot(v0,g0,v0)) == 0 # should be strictly real - @test abs(imag(dot(v0,g3,v0))) < 1e-10 - @test dot(v0,g0,v0) ≈ 0.65 rtol=0.01 - @test dot(v0,g1,v0) ≈ 0.67 rtol=0.01 - @test dot(v0,g2,v0) ≈ 0.67 rtol=0.01 - @test dot(v0,g3,v0) ≈ 0.65 rtol=0.01 - @test num_offdiagonals(g0,aIni) == 2 - - # on first component - g0f = G2MomCorrelator(0,:first) - g1f = G2MomCorrelator(1,:first) - @test imag(dot(v0,g0f,v0)) == 0 # should be strictly real - @test dot(v0,g0f,v0) ≈ 1.33 rtol=0.01 - @test dot(v0,g1f,v0) ≈ 1.33 + 7.08e-5im rtol=0.01 - # on second component - g0s = G2MomCorrelator(0,:second) - g1s = G2MomCorrelator(1,:second) - #@test_throws ErrorException("invalid ONR") get_offdiagonal(g0s,aIni,1) # should fail due to invalid ONR - @test dot(v0,g0s,v0) ≈ 1/3 - @test dot(v0,g1s,v0) ≈ 1/3 - # test against BoseFS - ham1 = HubbardMom1D(bfs1) - ham2 = HubbardMom1D(bfs2) - @test num_offdiagonals(g0f,aIni) == num_offdiagonals(ham1,bfs1) - @test num_offdiagonals(g0s,aIni) == num_offdiagonals(ham2,bfs2) - aIni = BoseFS2C(bfs2,bfs1) # flip bfs1 and bfs2 - @test get_offdiagonal(g0s,aIni,1) == (BoseFS2C(BoseFS{1,3}((0, 1, 0)),BoseFS{2,3}((1, 0, 1))), 0.47140452079103173) - # test on BoseFS - @test diagonal_element(g0s,bfs1) == 4/3 - @test diagonal_element(g0s,bfs2) == 1/3 - end - @testset "SuperfluidCorrelator" begin m = 6 n1 = 4 @@ -1231,7 +965,7 @@ using Rimu.Hamiltonians: circshift_dot @test diagonal_element(Momentum(1), BoseFS((1,0,0,0))) ≡ -1.0 @test_throws MethodError diagonal_element(Momentum(2), BoseFS((0,1,0))) - for address in (BoseFS2C((0,1,2,3,0), (1,2,3,4,5)), FermiFS2C((1,0,0,1), (0,0,1,0))) + for address in (FermiFS2C((1,0,0,1), (0,0,1,0)),) @test diagonal_element(Momentum(1), address) + diagonal_element(Momentum(2), address) ≡ diagonal_element(Momentum(0), address) end @@ -1247,7 +981,6 @@ using Rimu.Hamiltonians: circshift_dot for address in ( CompositeFS(BoseFS((1,2,3,4,5)), BoseFS((5,4,3,2,1))), - BoseFS2C((1,2,3,4,5), (5,4,3,2,1)) ) for i in 1:5 @test diagonal_element(DensityMatrixDiagonal(i, component=1), address) == i @@ -1574,7 +1307,6 @@ end @testset "TimeReversalSymmetry" begin @test_throws ArgumentError TimeReversalSymmetry(HubbardMom1D(BoseFS((1, 1)))) - @test_throws ArgumentError TimeReversalSymmetry(BoseHubbardMom1D2C(BoseFS2C((1, 1),(2,1)))) @test_throws ArgumentError begin TimeReversalSymmetry(HubbardRealSpace(CompositeFS(FermiFS((1, 1)),BoseFS((2,1))))) end @@ -1594,22 +1326,6 @@ end @test issymmetric(even_m) @test issymmetric(odd_m) end - @testset "2-particle BoseHubbardMom1D2C" begin - ham = BoseHubbardMom1D2C(BoseFS2C((0,1,1),(1,0,1))) - even = TimeReversalSymmetry(ham) - odd = TimeReversalSymmetry(ham; even=false) - - h_eigs = eigvals(Matrix(ham)) - p_eigs = sort!(vcat(eigvals(Matrix(even)), eigvals(Matrix(odd)))) - - @test starting_address(even) == time_reverse(starting_address(ham)) - @test h_eigs ≈ p_eigs - - @test issymmetric(Matrix(odd)) - @test issymmetric(Matrix(even)) - @test LOStructure(odd) isa IsHermitian - end - end @testset "Stoquastic" begin @@ -1849,7 +1565,7 @@ end @testset "dimension and multi-component addresses" begin addresses = [CompositeFS(FermiFS((1,0,1)), FermiFS((0,1,0))), BoseFS((1,0,1)), - FermiFS2C((1,0,1), (0,1,0)), BoseFS2C((1,0,1), (0,1,0)) + FermiFS2C((1,0,1), (0,1,0)) ] [@test dimension(addr) == dimension(typeof(addr)) for addr in addresses] end @@ -1928,7 +1644,7 @@ end t1 = 0; t2 = 0 for i in 1:4, j in i+1:4; t1 += 1; t2 = 0; - for k in 1:4, l in k+1:4 + for k in 1:4, l in k+1:4 t2+=1; tpd_f[t1,t2] = dot(dvec_f, TwoParticleExcitation(i, j, k, l), dvec_f); end end @@ -1939,5 +1655,8 @@ end @test LOStructure(op) isa IsHermitian test_observable_interface(ReducedDensityMatrix(1), BoseFS{4,4}(2,2,0,0)) test_observable_interface(ReducedDensityMatrix(2), FermiFS{2,4}(1,1,0,0)) + for r in (ReducedDensityMatrix(1), ReducedDensityMatrix{ComplexF32}(2)) + # Check that the result of show can be pasted into the REPL + @test eval(Meta.parse(repr(r))) == r + end end - diff --git a/test/allocation_tests.jl b/test/allocation_tests.jl index 902f793a3..026da2824 100644 --- a/test/allocation_tests.jl +++ b/test/allocation_tests.jl @@ -39,14 +39,6 @@ end ExtendedHubbardReal1D(b2), ExtendedHubbardReal1D(b3), - BoseHubbardReal1D2C(BoseFS2C(b1, b1)), - BoseHubbardReal1D2C(BoseFS2C(b2, b2)), - BoseHubbardReal1D2C(BoseFS2C(b3, b3)), - - BoseHubbardMom1D2C(BoseFS2C(b1, b1)), - BoseHubbardMom1D2C(BoseFS2C(b2, b2)), - BoseHubbardMom1D2C(BoseFS2C(b3, b3)), - HubbardRealSpace(b1), HubbardRealSpace(b2), HubbardRealSpace(b3), diff --git a/test/lomc.jl b/test/lomc.jl index 51eed21c1..ae1cb62dc 100644 --- a/test/lomc.jl +++ b/test/lomc.jl @@ -502,8 +502,7 @@ Random.seed!(1234) ) for address in ( - BoseFS2C((1,2,3), (0,1,0)), - CompositeFS(BoseFS((1,2,3)), FermiFS((0,1,0))) + CompositeFS(BoseFS((1,2,3)), FermiFS((0,1,0))), ) @test single_particle_density(address) == (1, 3, 3) @test single_particle_density(address; component=1) == (1, 2, 3) @@ -517,8 +516,7 @@ end @testset "Ground state energy estimates" begin for H in ( HubbardReal1D(BoseFS((1,1,2))), - BoseHubbardReal1D2C(BoseFS2C((1,2,2), (0,1,0))), - BoseHubbardMom1D2C(BoseFS2C((0,1), (1,0))), + HubbardMom1D(BoseFS((1, 1, 2))) ) @testset "$H" begin dv = DVec(starting_address(H) => 2; style=IsDynamicSemistochastic()) diff --git a/test/runtests.jl b/test/runtests.jl index dd9cf3038..a566bc46b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -97,51 +97,6 @@ using Rimu: replace_keys, delete_and_warn_if_present, clean_and_warn_if_others_p end end -@testset "BoseFS2C" begin - bfs2c = BoseFS2C(BoseFS((1,2,0,4)),BoseFS((4,0,3,1))) - @test typeof(bfs2c) <: BoseFS2C{7,8,4} - @test num_occupied_modes(bfs2c.bsa) == 3 - @test num_occupied_modes(bfs2c.bsb) == 3 - @test onr(bfs2c.bsa) == [1,2,0,4] - @test onr(bfs2c.bsb) == [4,0,3,1] - @test Hamiltonians.bose_hubbard_2c_interaction(bfs2c) == 8 # n_a*n_b over all sites -end - -@testset "TwoComponentBosonicHamiltonian" begin - aIni2cReal = BoseFS2C(BoseFS((1,1,1,1)),BoseFS((1,1,1,1))) # real space two-component - Ĥ2cReal = BoseHubbardReal1D2C(aIni2cReal; ua = 6.0, ub = 6.0, ta = 1.0, tb = 1.0, v= 6.0) - hamA = HubbardReal1D(BoseFS((1,1,1,1)); u=6.0, t=1.0) - hamB = HubbardReal1D(BoseFS((1,1,1,1)); u=6.0) - @test hamA == Ĥ2cReal.ha - @test hamB == Ĥ2cReal.hb - @test num_offdiagonals(Ĥ2cReal,aIni2cReal) == 16 - @test num_offdiagonals(Ĥ2cReal,aIni2cReal) == num_offdiagonals(Ĥ2cReal.ha,aIni2cReal.bsa)+num_offdiagonals(Ĥ2cReal.hb,aIni2cReal.bsb) - @test dimension(Ĥ2cReal) == 1225 - @test dimension(Float64, Ĥ2cReal) == 1225.0 - - hp2c = offdiagonals(Ĥ2cReal,aIni2cReal) - @test length(hp2c) == 16 - @test hp2c[1][1] == BoseFS2C(BoseFS((0,2,1,1)), BoseFS((1,1,1,1))) - @test hp2c[1][2] ≈ -1.4142135623730951 - @test diagonal_element(Ĥ2cReal,aIni2cReal) ≈ 24.0 # from the V term - - aIni2cMom = BoseFS2C(BoseFS((0,4,0,0)),BoseFS((0,4,0,0))) # momentum space two-component - Ĥ2cMom = BoseHubbardMom1D2C(aIni2cMom; ua = 6.0, ub = 6.0, ta = 1.0, tb = 1.0, v= 6.0) - @test num_offdiagonals(Ĥ2cMom,aIni2cMom) == 9 - @test dimension(Ĥ2cMom) == 1225 - @test dimension(Float64, Ĥ2cMom) == 1225.0 - - hp2cMom = offdiagonals(Ĥ2cMom,aIni2cMom) - @test length(hp2cMom) == 9 - @test hp2cMom[1][1] == BoseFS2C(BoseFS((1,2,1,0)), BoseFS((0,4,0,0))) - @test hp2cMom[1][2] ≈ 2.598076211353316 - - smat2cReal, adds2cReal = ExactDiagonalization.build_sparse_matrix_from_LO(Ĥ2cReal,aIni2cReal) - eig2cReal = eigen(Matrix(smat2cReal)) - smat2cMom, adds2cMom = ExactDiagonalization.build_sparse_matrix_from_LO(Ĥ2cMom, aIni2cMom) - eig2cMom = eigen(Matrix(smat2cMom)) - @test eig2cReal.values[1] ≈ eig2cMom.values[1] -end @safetestset "KrylovKit" begin include("KrylovKit.jl")