Add FermiFCI-1.0

FermiFCI.jl-3a4bc10261aa112194195c82333a9f6ed312c004/.gitignore

+dev_scripts/
+*output/
+*_archive/

FermiFCI.jl-3a4bc10261aa112194195c82333a9f6ed312c004/Compat.toml

+[0]
+julia="1.5.3-*"

FermiFCI.jl-3a4bc10261aa112194195c82333a9f6ed312c004/LICENSE

+Copyright 2021 Lukas Rammelmüller, Artem G. Volosniev, David Huber
+
+Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

FermiFCI.jl-3a4bc10261aa112194195c82333a9f6ed312c004/Project.toml

+name = "FermiFCI"
+uuid = "c98d26ac-406e-49de-a033-695b2b85f5d1"
+authors = ["ramelmueller <lukas.rammelmueller@gmail.com>"]
+version = "0.1.0"
+
+[deps]
+Arpack = "7d9fca2a-8960-54d3-9f78-7d1dccf2cb97"
+Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
+DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
+DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
+HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
+LoggingExtras = "e6f89c97-d47a-5376-807f-9c37f3926c36"
+PProf = "e4faabce-9ead-11e9-39d9-4379958e3056"
+PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0"
+QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
+Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
+SpecialPolynomials = "a25cea48-d430-424a-8ee7-0d3ad3742e9e"
+YAML = "ddb6d928-2868-570f-bddf-ab3f9cf99eb6"

FermiFCI.jl-3a4bc10261aa112194195c82333a9f6ed312c004/README.md

+# FermiFCI.jl
+A simple and lightweight toolbox for performing exploratory FCI calculation in arbitrary single-particle bases written in the Julia language.
+
+## Getting started
+Here's a few tips on how to get started, assuming you have downloaded Julia already. If you haven't please check out [the documentation on how to get started with Julia](https://docs.julialang.org/en/v1/manual/getting-started/).
+
+### Prerequisites
+The package was developed with `Julia v1.5.3` so make sure you have this or a later version installed on your machine.
+
+### Installing `FermiFCI.jl`
+To get started with `FermiFCI.jl` simply clone this repository and add the package to your global or local scope. Generally, it is recommended to use a local scope for each of your Julia projects. This can be done by starting the Julia REPL as follows:
+```
+    julia --project=/path/to/the/root/of/your/project
+```
+When started this way, newly installed packages are added only to this project. This is particularly helpful on HPC environments as it circumvents potential permission issues. Packages can be added either directly in the Julia REPL, or by entering the Pkg REPL with `]` and simply type `add PackageName` or, for the case here, `add /path/to/FermiFCI/`. More information can be found [here](https://docs.julialang.org/en/v1/stdlib/Pkg/).
+
+Once you added the `FermiFCI.jl` package, this is it - you're good to go.
+
+Roadmap: In the near future FermiFCI will be available as an official package - stay tuned!
+
+
+### Usage
+The package is deliberately kept simple - there's only a handful of routines to call, mainly to construct a hamiltonian and to diagonalize the resulting matrix along with some type definitions and a few extra utility functions to carry out common tasks (such as the construction of a plain list of available states).
+
+As a first step, it is advisable to execute one of the examples which are discussed [in this paper]() - the paper also explains some generic FCI details. The source code of those is available in a [separate repository](https://github.com/rammelmueller/fermifci_data_repo). Here, it's good to know that adding the package itself only installs the required prerequisites for the core functionality of `FermiFCI.jl`. Packages that are used by the examples need to be added manually - see the list of required packages in the documentation of the examples themselves.
+
+Should you choose not too go with one of the provided examples (or moved on from this first step), all you have to do is to include `using FermiFCI` in your code, then you'll be able to use all core functionality.
+
+
+### Troubleshooting
+- It was brought to our attention that on MacOS some imports require adjusting filepaths to absolute paths, which seems to be connected with the HDF5.jl package.
+
+
+## Applications
+Here we list some applications that rely on FermiFCI.jl:
+- [A modular implementation of an effective interaction approach for harmonically trapped fermions in 1D](https://arxiv.org/abs/2202.04603) - accompanying paper to this package. Implements harmonically trapped fermions with and without energy restriction of the many-body basis, few-fermion problems with flat traps with and without mass-imbalance of the two species as well as an effective interaction approach.
+- [Magnetic impurity in a one-dimensional few-fermion system](https://arxiv.org/abs/2204.01606) - Extended treatment of harmonically trapped fermions with the presence of a delta-like static magnetic impurity.
+
+If your paper should be included in this list, please get in touch!
+
+## Legal
+Published under the MIT Licencse. If you use this for any publication, please cite us. Have fun!

FermiFCI.jl-3a4bc10261aa112194195c82333a9f6ed312c004/src/FermiFCI.jl

+#===============================================================================
+
+    FermiFCI.jl - LR, June 2021
+
+===============================================================================#
+module FermiFCI
+
+# Those are the exported functions.
+# (the ones that can be called without the namespace prefix)
+export Orbital, FullState, SpinState, OneBodyCoeffTensor, TwoBodyCoeffTensor, construct_hamiltonian, diagonalize, mb_state_energy
+
+
+
+# Holds type definitions, should be loaded first.
+include("typedefs.jl")
+include("state_reps/state_functions.jl")
+
+# Some utilities for basis preparation.
+include("utils/plain_hilbert_construction.jl")
+
+# Density matrices.
+include("utils/density_matrix.jl")
+include("utils/density_profile.jl")
+
+# Stuff for constructing the actual Hamiltonian.
+include("construction.jl")
+
+# Stuff for diagonalization.
+include("diagonalization.jl")
+
+# Define a submodule IO.
+module Utils
+    include("io/io_helper.jl")
+end # module IO
+
+end # module FermiFCI

FermiFCI.jl-3a4bc10261aa112194195c82333a9f6ed312c004/src/README.md

+## Core functionality
+The core functionality itself should not be altered by the user. This explanation simply provides some insights into the inner workings of the package and which datatypes are available.
+
+
+### Data types
+`FermiFCI.jl` provides a few data types to represent physical objects. These are mostly self-explanatory and (mostly) defined (as well as documented) in the file `src/typedefs.jl`. Note, that most of these data types are merely type aliases to provide better readablity as well as a more abstract implementation.
+
+
+### Representation of states
+The physical states are represented as binary numbers where each bit represents the occupation of a given single-particle orbital (which, due to the fermions, is either 0 or 1). For instance, if we want to represent a many-body state with particles of a given spin in the first and fourth single-particle orbitals, we may do so with the bitstring `01001` (to be read from right to left). To accomodate this, `FermiFCI.jl` provides the data type `SpinState` - all objects that represent the many-body state of a single spin species should be of this type. The previous example state, for instance, may be produced by first constructing an empty state and then creating particles in the corresponding orbitals:
+```
+    state = SpinState(0) # Empty state.
+    state = create(state, 1) # Add first particle.
+    state = create(state, 4) # Add second particle.
+```
+which uses the provided functions for creation operators (there's a corresponding annihilation function). Alternatively, we could simply cast from the integer value: `state = SpinState(9)`.
+
+While the current implementation uses simple unsigned integer types, this may change in the future. Moreover, the abstract nature of the implementation allows to quickly switch to an alternate (and custom) representation, if required.
+
+To reflect a full state with both spin species, two objects of the `SpinState` type are combined to a `FullState`. To this end, the conversion routine `f = s_to_f(s, m)` is provided which converts two objects of the type `SpinState` to a `FullState` object. Accordingly, the expression `s, m = f_to_s(f)` does the opposite.
+
+In addition, these states have some extra functions defined which may be looked up directly in the source code in `src/state_reps/`.
+
+
+### Construction of the Hamiltonian
+The function `construct_hamiltonian` is the main tool of this package. It's function is to simply compute all matrix elements of the Hamiltonian of interest and return an object of the type `Hamiltonian`, which essentially represents a sparse matrix representation of the full Hamiltonian matrix (and therefore requires some memory for large systems). The function is generic for all two-species Hamiltonians - the physical content is encoded in the provided coefficient tensors.
+
+To use this function to obtain all matrix elements, one may simply use 
+```
+    ham = construct_hamiltonian(mb_basis; up_coeffs, down_coeffs, up_down_coeffs, up_up_coeffs, down_down_coeffs)
+```
+where the argument `mb_basis` reflects the many-body hilbert space which is nothing but a list of states of type `FullState`. The remaining arguments are optional and represent the corresponding single- and two-body coefficients of the type `OneBodyCoeffTensor` and `TwoBodyCoeffTensor`, respectively.
+
+
+### Diagonalization
+Once the elements of a given Hamiltonian have been computed, all that is left to be done is the diagonalization of the matrix. This can be achieved via calling the function `diagonalize` - this is essentially a wrapper for the ARPACK function `eigs` which uses the iterated Arnoldi method. For the previously constructed Hamiltonian `ham`, for example, this would look like
+```
+    ev, est = diagonalize(ham, param)
+```
+where extra parameters such as the number of eigenvalues (`n_eivenvalues`), the type of eigenvalues (`ev_type`), the absolute tolerance (`tol`) as well as the maximum number of iterations (`max_iter`) may be specified in the dictionary `param` with the indicated keys. If nothing is specified default values are used. All these parameters are described [in the ARPACK.jl package](https://arpack.julialinearalgebra.org/latest/).
+
+The diagonalization returns the final result of the calculation: Eigenvalues (`ev`) and eigenstates (`est`) for the constructed Hamiltonian. 
+
+
+## Utils
+These are various functions that are not strictly speaking in the core functionality of FermiFCI. Examples include the comuptation of the one-body density matrix, density profiles and construction functions for the many body basis. Those are found in `utils/`.

FermiFCI.jl-3a4bc10261aa112194195c82333a9f6ed312c004/src/construction.jl

+#===============================================================================
+
+    construction.jl - LR, June 2020
+
+    Routines to actually construct the Hamiltonian from a given many-body basis.
+
+===============================================================================#
+using SparseArrays, LinearAlgebra
+include("typedefs.jl")
+
+
+function get_particle_hops(state::SpinState, n_basis::Integer)::Array{Tuple{Int,Int,SpinState,Int},1}
+    """ Returns possible single particle jumps from one state to another. This
+        requires relatively little memory but using this speeds up the computation
+        slightly and we don't have to re-compute things later (for additional
+        terms, for instance).
+
+        Number of entries per state: n_part*(n_basis-n_npart+1)
+        Example: n_basis = 30 // n_part = 16 ---> 240 single particle hops
+
+    """
+    hops = Array{Tuple{Int,Int,SpinState,Int},1}()
+    for i = 1:n_basis
+        temp = annihilate(state, i)
+        if !isnothing(temp)
+            for j = 1:n_basis
+                new_state = create(temp, j)
+                if !isnothing(new_state)
+                    push!(hops, (i,j,new_state,sign(state,i,j)))
+                end
+            end
+        end
+    end
+    return hops
+end
+
+
+function find_n_basis(hilbert_space::Array{FullState,1})::Integer
+    """ Finds the maximal index of the occupied single-partilce orbitals so that
+        the loops are terminated at this level.
+
+        Example:
+
+                states = [|1000000>,|0000100>] -> find_n_basis(states) = 7
+
+        Notes:
+            -   Potentially this strategy could be produce a hit on performance,
+                but it is also the most general one.
+            -   TODO: this could be optimized if the lookup_table is assumed to
+                be ordered - then the last entry would automatically give the
+                correct answer -> would need OrderedDict. Is this optimization
+                necessary?
+    """
+    n_basis::Integer = 0
+    for state in hilbert_space
+        s_up, s_down = f_to_s(state)
+
+        for s in [s_up, s_down]
+            m = max_orbital(s)
+            if m > n_basis
+                n_basis = m
+            end
+        end
+    end
+    return n_basis
+end
+
+
+function construct_hamiltonian(
+        hilbert_space::Array{FullState,1};
+        up_coeffs::Union{Nothing,OneBodyCoeffTensor}=nothing, # ↑ single-body coefficients
+        down_coeffs::Union{Nothing,OneBodyCoeffTensor}=nothing, # ↓ single-body coefficients
+        up_down_coeffs::Union{Nothing,TwoBodyCoeffTensor}=nothing, # ↑↓ interaction
+        up_up_coeffs::Union{Nothing,TwoBodyCoeffTensor}=nothing, # ↑↑ interaction
+        down_down_coeffs::Union{Nothing,TwoBodyCoeffTensor}=nothing # ↓↓ interaction
+    )::Hamiltonian where {T <: Orbital}
+    """ Loops through the entire Hilbert space and finds the Hamiltonian Matrix
+        elements. Returns a Hamiltonian structure.
+    """
+    # These values are the entries of the Hamiltonian. Assuming that they wont
+    # be larger than 2^32, we can use 32 bit integers here.
+    row = Vector{SparseIndexType}()
+    col = Vector{SparseIndexType}()
+    data = Vector{DType}()
+
+    # Find maximally occupied index.
+    n_basis::Integer = find_n_basis(hilbert_space)
+    lookup_table, inv_lookup_table = make_lookup_table(hilbert_space)
+
+    # Loop over all states in the Hilbert space.
+    for n = 1:length(lookup_table)
+        # Split into single spin states.
+        s_up, s_down = f_to_s(lookup_table[n])
+
+        up_hops = get_particle_hops(s_up, n_basis)
+        down_hops = get_particle_hops(s_down, n_basis)
+
+        # ----------------------
+        # One-body terms.
+
+        if !isnothing(up_coeffs)
+            for (l,i,new_up,sign) in up_hops
+                new_state = get(inv_lookup_table, s_to_f(new_up, s_down), nothing)
+                if !isnothing(new_state)
+                    push!(col, n)
+                    push!(row, new_state)
+                    push!(data, -sign * up_coeffs[i,l])
+                end
+            end
+        end # End  of ↑ section.
+
+        if !isnothing(down_coeffs)
+            for (l,i,new_down,sign) in down_hops
+                new_state = get(inv_lookup_table, s_to_f(s_up, new_down), nothing)
+                if !isnothing(new_state)
+                    push!(col, n)
+                    push!(row, new_state)
+                    push!(data, -sign * down_coeffs[i,l])
+                end
+            end
+        end # End  of ↓ section.
+
+        # ----------------------
+
+        # Two-body part.
+        # Loop structure:
+        #   k -> mu anihilator
+        #   i -> mu creator
+        #   l -> nu anihilator
+        #   j -> nu creator
+
+        # UP/DOWN loop (mu =  ↑, nu = ↓).
+        if !isnothing(up_down_coeffs)
+            for (k,i,new_up,sign_up) in up_hops
+                for (l,j,new_down,sign_down) in down_hops
+                    new_state = get(inv_lookup_table, s_to_f(new_up, new_down), nothing)
+                    if !isnothing(new_state)
+                        push!(col, n)
+                        push!(row, new_state)
+                        push!(data, sign_up*sign_down * up_down_coeffs[i,j,k,l])
+                    end
+                end
+            end
+        end # End  of ↑↓ section.
+
+        # UP/UP loop  (mu =  ↑, nu = ↑).
+        if !isnothing(up_up_coeffs)
+            for (k,i,temp_up,temp_sign) in up_hops
+                temp_hops = get_particle_hops(temp_up, n_basis)
+                for (l,j,new_up,new_sign) in temp_hops
+                    new_state = get(inv_lookup_table, s_to_f(new_up, s_down), nothing)
+                    if !isnothing(new_state)
+                        push!(col, n)
+                        push!(row, new_state)
+                        push!(data, temp_sign*new_sign * up_up_coeffs[i,j,k,l])
+                    end
+                end
+            end
+        end # End  of ↑↑ section.
+
+        # DOWN/DOWN loop  (mu =  ↓, nu = ↓).
+        if !isnothing(down_down_coeffs)
+            for (k,i,temp_down,temp_sign) in down_hops
+                temp_hops = get_particle_hops(temp_down, n_basis)
+                for (l,j,new_down,new_sign) in temp_hops
+                    new_state = get(inv_lookup_table, s_to_f(s_up, new_down), nothing)
+                    if !isnothing(new_state)
+                        push!(col, n)
+                        push!(row, new_state)
+                        push!(data, temp_sign*new_sign * down_down_coeffs[i,j,k,l])
+                    end
+                end
+            end
+        end # End  of ↓↓ section.
+
+
+        # ----------------------
+
+    end # End of state loop.
+
+    # Return a Hamiltonian object.
+    return Hamiltonian(row, col, data, length(lookup_table))
+end

FermiFCI.jl-3a4bc10261aa112194195c82333a9f6ed312c004/src/diagonalization.jl

+using SparseArrays, LinearAlgebra, Arpack
+include("io/io_helper.jl")
+
+
+function diagonalize(ham::Hamiltonian, param::Dict{Any,Any})
+    """ Performs the diagonalization as specified in the paramter dictionary.
+        Done either by iterated Arnoldi method with ARPACK (when full_diag is
+        set to `false`) or by full diagonalization.
+
+        Notes:
+         -  Julia has a matrix-independent iteration number of 300 per default,
+            this is potentially problematic. Needs to be adjusted in some cases.
+    """
+    sparse_ham = sparse(ham.row, ham.col, ham.data, ham.n_fock, ham.n_fock)
+
+    if ham.n_fock > 1
+        @info "Starting to diagonalize the Hamiltonian." n_fock=Int(ham.n_fock)
+        if get(param, "full_diag", false)
+            # Full diagonalization with built-in method.
+            time = @elapsed mem = @allocated F = eigen(Matrix(sparse_ham))
+            @info "Done computing the full spectrum." time=time allocated=MemoryTag(mem) spectrum=F.values
+            return F.values, F.vectors
+        else
+            # Arnoldi method with ARPACK.
+            which_map = Dict([
+                ("SA", :SR)
+            ])
+            # Here, we set the default values that are otherwise overwritten
+            # if specified in the parameter dictionary.
+            time = @elapsed mem = @allocated ev, est = eigs(
+                sparse_ham,
+                nev=get(param, "n_eigenvalues", 10),
+                which=which_map[get(param, "ev_type", "SA")],
+                tol=get(param, "tol", 1e-15),
+                maxiter=get(param, "max_iter", 10*ham.n_fock)
+            )
+            @info "Done computing the lower spectrum." time=time allocated=MemoryTag(mem) spectrum=ev
+            return ev, est
+        end
+    end
+    @info "Trivial 1x1 matrix to diagonalize."
+    return Matrix(sparse_ham)[1,1], [1]
+end

FermiFCI.jl-3a4bc10261aa112194195c82333a9f6ed312c004/src/io/io_helper.jl

+#===============================================================================
+
+    io_helper.jl - LR, June 2021
+
+    Some useful functionality for nicer output.
+
+===============================================================================#
+struct MemoryTag
+    """ Human readable memory tag.
+    """
+    bytes::Integer
+end
+
+function Base.show(io::IO, tag::MemoryTag)
+    print(io, round(tag.bytes/1024^2, digits=3), "MB")
+end

FermiFCI.jl-3a4bc10261aa112194195c82333a9f6ed312c004/src/state_reps/rep_integer.jl

+#===============================================================================
+
+    rep_integer.jl - LR, Jan 2021
+
+    Representation of states is done via an Integer.
+
+===============================================================================#
+const SpinState = UInt64 # That's a single species state.
+const FullState = UInt128 # That's the full state describing the entire system.
+
+# That's the shift which is applied to fuse the two spins to a full state.
+const _bshift = 2^60
+
+function f_to_s(fs::FullState)::Tuple{SpinState,SpinState}
+    """ Takes a state that represents the full state in the Hilbert space
+        and maps it to the constituent spin states.
+    """
+    return SpinState(fs%_bshift), SpinState(fs÷_bshift)
+end
+
+function s_to_f(su::SpinState, sd::SpinState)::FullState
+    """ Maps two spin states to a full state.
+    """
+    return FullState(sd)*_bshift + FullState(su)
+end
+
+
+function create(s::SpinState, i::Integer)::Union{Nothing,SpinState}
+    """ Creates a particle at position i, false if there's a particle already.
+    """
+    shift = 1 << (i-1)
+    if s & shift == 0
+        return s ⊻ shift
+    else
+        return nothing
+    end
+end
+
+function annihilate(s::SpinState, i::Integer)::Union{Nothing,SpinState}
+    """ Destroys a particle at position i, false if there's no particle to destroy.
+    """
+    shift = 1 << (i-1)
+    if s & shift > 0
+        return s ⊻ shift
+    else
+        return nothing
+    end
+end
+
+function count_occupancies(s::SpinState, src::Integer, dst::Integer)
+    # return sum([(s >>> k) & 1 for k=0:63]) --> slower!
+    c = 0
+    for k=src-1:dst-1
+        c += (s >>> k) & 1
+    end
+    return c
+end
+
+function sign(s::SpinState, src::Integer, dst::Integer)
+    """ This is the sum of densities between the anihilation and creation operators
+        (including both endpoints) which gives the entire fermionic sign factor.
+        The contributions from up and down states simply add.
+    """
+    if  src > dst
+        return (-1)^count_occupancies(s, dst, src)
+    end
+    return (-1)^count_occupancies(s, src, dst)
+end
+
+
+
+# Indexing.
+function max_orbital(s::SpinState)
+    """ Returns the index of the maximally occupied orbital in a SpinState.
+
+        Note:
+            - For integers, this is essentially the bit length.
+    """
+    # sizeof = number of bytes.
+    # leading_zeros = built in for leading zeros in representation.
+    return sizeof(s) * 8 - leading_zeros(s)
+end
+function Base.getindex(s::SpinState, i::Unsigned)::SpinState
+    """ Note: Indexing works with the lowest bits first (i.e., the ones from the
+        right). It doesn't matter how it is used, only the implementation here
+        needs to be consistent.
+
+        Only positive integers are allowed. Negatives would only give zeros.
+    """
+    return s >>> (i-1) & 1
+end
+Base.getindex(s::SpinState, i::Integer)::SpinState = s[UInt(i)]
+Base.firstindex(s::SpinState) =  1
+Base.lastindex(s::SpinState) = max_orbital(s)
+Base.length(s::SpinState) = max_orbital(s)

FermiFCI.jl-3a4bc10261aa112194195c82333a9f6ed312c004/src/state_reps/state_functions.jl

+#===============================================================================
+
+    state_functions.jl
+
+    Some functions that states support regardless of the specific implementation.
+
+===============================================================================#
+function mb_state_energy(orbital::T, state::SpinState)::AbstractFloat where T<:Orbital
+    """ Retrieves the energy for a many-body state (single spin species).
+    """
+    e = 0
+    for k = 1:length(state)
+        if state[k]>0
+            e += orbital(k)
+        end
+    end
+    return e
+end
+
+
+function mb_state_energy(orbital_up::T, orbital_down::S, state::FullState)::AbstractFloat where {T<:Orbital,S<:Orbital}
+    """ Many-body energy for a full state (wraps the single-particle routine)..
+    """
+    state_up, state_down = f_to_s(state)
+    return mb_state_energy(orbital_up, state_up) + mb_state_energy(orbital_down, state_down)
+end

FermiFCI.jl-3a4bc10261aa112194195c82333a9f6ed312c004/src/typedefs.jl

+#===============================================================================
+
+    typedefs.jl - LR, June 2020
+
+    Holds some definitions of types.
+
+===============================================================================#
+const DType = Float64
+const CType = Complex{DType}
+const IType = Int64
+
+# Types for state handling. The type of representation is chosen by including the
+# corresponding file.
+include("state_reps/rep_integer.jl")
+
+# Index in the single-particle basis
+# Notes:
+#   -   Limited to 256 for now, can be extended arbitrarily.
+#   -   Don't make unsigned: differenes (which are needed!) will also be unsigned
+#       then -> this does prevent us from checking the positivity.
+# const BasisIndex = Int16
+
+# Index in the many-body Hilbert space.
+# Notes:
+#   -   This is the dimension of the many-body Hilbert space, a few millions
+#       should be reachable.
+const HilbertIndex = UInt32
+
+# Default datatypes for the lookup dictionaries.
+const LookupDict = Dict{HilbertIndex,FullState}
+const InvLookupDict = Dict{FullState,HilbertIndex}
+
+# Type for index in sparse matrices.
+# Notes:
+#   -   This must be at least as large as HilbertIndex, however, since there will
+#       be many entries UInt64 is recommended for large runs (with more than dim H = 1e6).
+#   -   This requirement is a bit unintuitive at first, but it has to do how the
+#       sparse matrices are stored (compressed format, see the CSC standard).
+#   -   Naturally, a larger datatype will increase memory consumption, but this is
+#       what we have to live with for larger runs.
+const SparseIndexType = UInt64
+
+
+# Types for storage of numerical coefficients.
+const WaveFunction{T<:Union{DType,CType}} = Array{T,1}
+const OneBodyCoeffTensor{T<:DType} = Array{T,2}
+const TwoBodyCoeffTensor{T<:DType} = Array{T,4}
+
+abstract type Orbital end
+
+
+struct Hamiltonian
+    """ Datatype that represents the Hamiltonian. Merely holds all necessary
+        information to construct a sparse matrix and is used only for better
+        readability.
+    """
+    row::Vector{HilbertIndex}
+    col::Vector{HilbertIndex}
+    data::Vector{DType}
+    n_fock::HilbertIndex
+end

FermiFCI.jl-3a4bc10261aa112194195c82333a9f6ed312c004/src/utils/density_matrix.jl

+include("../construction.jl")
+
+
+function compute_obdm(wf::WaveFunction, flavor::Integer, hilbert_space::Array{FullState,1})::Array{DType,2}
+    """ Constructs the one-body density matrix for a given wavefunction.
+    """
+    n_basis::Integer = find_n_basis(hilbert_space)
+    lookup_table, inv_lookup_table = make_lookup_table(hilbert_space)
+
+    obdm = zeros(DType, (n_basis, n_basis))
+    for n = 1:length(lookup_table)
+        state = lookup_table[n]
+
+        # Choose up or down component.
+        (state, rest) = f_to_s(state)[[flavor,3-flavor]]
+        for k = 1:n_basis
+            temp = annihilate(state, k)
+            if !isnothing(temp)
+                for i = 1:n_basis
+                    new_state = create(temp, i)
+                    if !isnothing(new_state)
+                        full_state = (flavor == 1) ? s_to_f(new_state, rest) : s_to_f(rest, new_state)
+
+                        # This is the index of the final state.
+                        final_state = get(inv_lookup_table, full_state, nothing)
+                        if !isnothing(final_state)
+                            obdm[k,i] += sign(state, k, i) * wf[n]*conj(wf[final_state])
+                        end
+                    end
+                end
+            end
+        end
+    end
+
+    return obdm
+end
+
+
+function compute_tbdm_up_down(wf::WaveFunction, hilbert_space::Array{FullState,1})::Array{DType,4}
+    n_basis::Integer = find_n_basis(hilbert_space)
+    lookup_table, inv_lookup_table = make_lookup_table(hilbert_space)
+
+    # The one-body density-matrix in the
+    tbdm = zeros(DType, (n_basis, n_basis, n_basis, n_basis))
+    for n = 1:length(lookup_table)
+        s_up, s_down = f_to_s(lookup_table[n])
+        for k = 1:n_basis
+            temp_up = annihilate(s_up, k)
+            if !isnothing(temp_up)
+                for i = 1:n_basis
+                    new_up = create(temp_up, i)
+                    if !isnothing(new_up)
+                        for l = 1:n_basis
+                            temp_down = annihilate(s_down, l)
+                            if !isnothing(temp_down)
+                                for j = 1:n_basis
+                                    new_down = create(temp_down, j)
+                                    if !isnothing(new_down)
+                                        final_state = inv_lookup_table[s_to_f(new_up,new_down)]
+                                        tbdm[k,i,l,j] += sign(s_up,k,i)*sign(s_down,l,j) * wf[n]*conj(wf[final_state])
+                                    end
+                                end
+                            end
+                        end
+                    end
+                end
+            end
+        end
+    end
+
+    return tbdm
+end

FermiFCI.jl-3a4bc10261aa112194195c82333a9f6ed312c004/src/utils/density_profile.jl

+function compute_density_profile(orbital::T, xgrid::Array{DType,1}, obdm::Array{DType,2}; kwargs...)::Array{DType,1} where T<:Orbital
+    """ Takes the one-body density matrix and returns the spatial density-profile
+        on the specified position grid (unnormalized).
+    """
+    density = zeros(DType, (length(xgrid),))
+    for i=1:size(obdm)[1]
+        di = conj.(orbital.(i, xgrid))
+        for j=1:size(obdm)[2]
+            density .+= di.* orbital.(j, xgrid) .* obdm[i,j]
+        end
+    end
+    return -density
+end
+
+
+function compute_updown_noise_correlation(orbital::T, xgrid::Array{DType,1}, tbdm::Array{DType,4}; kwargs...)::Array{DType,2} where T<:Orbital
+    """ Takes the two-body density matrix (for up and down) and returns the spatial noise-correlation-profile
+        on the specified position grid (unnormalized).
+    """
+    nn_corr = zeros(DType, (length(xgrid),length(xgrid)))
+    for i=1:size(obdm)[1]
+        di = conj.(orbital.(i, xgrid))
+        for j=1:size(obdm)[2]
+            density .+= di.* orbital.(j, xgrid) .* obdm[i,j]
+        end
+    end
+    return -density
+end

FermiFCI.jl-3a4bc10261aa112194195c82333a9f6ed312c004/src/utils/plain_hilbert_construction.jl

+#===============================================================================
+
+    plain_hilbert_construction.jl - LR, June 2020
+
+    Functionality for creating a plain cutoff basis, i.e., all many-body states
+    of given particle number within a single-particle basis up to some cutoff.
+    This is the most generic (however, mostly used) case.
+
+===============================================================================#
+
+
+function _find_states(state::SpinState, nb::Integer, n::Integer, l::Integer)::Array{SpinState,1}
+    """ Recursive piece of the state-finding algorithm. Terminates once the
+        maximum length is reached.
+    """
+    if n == 0
+        return [state << k for k = 0:(nb-l)]
+    end
+    if (nb - l) < n
+        return []
+    end
+    if l > 0
+        return [_find_states((state<<1)+1, nb, n-1, l+1); _find_states(state<<1, nb, n, l+1)]
+    end
+    return _find_states((state<<1)+1, nb, n-1, l+1)
+end
+
+# ------------------------------------------------------------------------------
+
+function get_plain_fock_basis(n_basis::Integer, n_part::Array{IType,1})::Array{FullState,1}
+    """ Constructs the many-body basis with a plain basis-state cutoff.
+
+        Attention: Currently limited to 64 bit.
+    """
+    s_up = _find_states(SpinState(0), n_basis, n_part[1], 0)
+    s_down = _find_states(SpinState(0), n_basis, n_part[2], 0)
+    return sort([s_to_f(u,d) for u in s_up for d in s_down])
+end
+
+
+# ------------------------------------------------------------------------------
+
+
+function make_lookup_table(states::Array{FullState,1})::Tuple{LookupDict,InvLookupDict}
+    """ Takes a list of states and produces the lookup tables.
+    """
+    lookup_table = LookupDict()
+    inverse_lookup_table = InvLookupDict()
+    for k = 1:length(states)
+        lookup_table[k] = states[k]
+        inverse_lookup_table[states[k]] = k
+    end
+
+    return lookup_table, inverse_lookup_table
+end