cc_amplitudes.F90 Source File


Contents

Source Code


Source Code

#include "macros.h"

module cc_amplitudes

    use SystemData, only: nbasis, nel, nOccAlpha, nOccBeta, ElecPairs
    use detbitops, only: get_bit_excitmat, FindBitExcitLevel
    use FciMCData, only: totwalkers, ilutref, currentdets, AllNoatHf, projedet, &
                         HashIndex, CurrentDets, ll_node
    use constants, only: dp, lenof_sign, EPS, n_int, bits_n_int, stdout
    use bit_rep_data, only: niftot, nifd, extract_sign
    use bit_reps, only: decode_bit_det
    use replica_data, only: AllEXLEVEL_WNorm
    use back_spawn, only: setup_virtual_mask, mask_virt_ni
    use hash, only: hash_table_lookup, FindWalkerHash
    use util_mod, only: swap, choose_i64, operator(.div.)
    use bit_rep_data, only: nifd
    use MPI_wrapper, only: iProcIndex, root
    use Parallel_neci, only: MPISumAll, MPIReduce, MPI_SUM, MPI_LOR, MPIAllLorLogical
    use util_mod, only: stop_all

    implicit none
    private
    public :: cc_amplitude, calc_number_of_excitations, &
        calc_n_parallel_excitations, setup_ind_matrix_singles, &
        t_plot_cc_amplitudes, t_cc_amplitudes, &
        cc_singles_factor, cc_doubles_factor, cc_triples_factor, cc_quads_factor, &
        cc_delay, init_cc_amplitudes, print_cc_amplitudes, cc_order

    logical :: t_cc_amplitudes = .false.
    logical :: t_plot_cc_amplitudes = .false.

    ! we need the order of cluster operators..
    integer :: cc_order = 0

    integer :: cc_delay = 1000

    integer :: est_triples = 0
    integer :: n_singles = 0, n_doubles = 0
    integer :: n_triples = 0, n_quads = 0
    ! also store the L0, L1 and L2 norm of the cc-amps up to quadrupels
    ! so i dont need the est_ quantities from above
    real(dp) :: cc_amp_norm(0:2, 4) = 0.0_dp

    ! and also create it for the mpi sum..
    integer :: all_est_triples = 0, all_n_singles = 0, all_n_doubles = 0
    real(dp) :: all_cc_amp_norm(0:2, 4) = 0.0_dp

    logical :: t_store_hash_quadrupels = .true.
    integer :: quad_hash_size, quad_hash_size_wf

    ! hard code for now which norm should be compared..
    integer :: norm_comp = 2
    ! i guess i need a new hash table for cc amplitudes, since i also
    ! want the amplitudes.. but how do i deal with clashes here?
    type cc_hash
        logical :: found = .false.
        real(dp) :: amp = 0.0_dp
        integer(n_int), allocatable :: ind(:)
        type(cc_hash), pointer :: next => null()
    end type cc_hash

    type(cc_hash), pointer :: quad_hash(:)

    type(cc_hash), pointer :: quad_hash_wf(:)

    integer :: n_clashes = 0

    ! maybe it would be nice to have a type which encodes this information
    ! and which gives an easy and nice way to de/encode the indices
    ! involved..
    type cc_amplitude
        integer :: order = 0
        integer, allocatable :: operators(:, :, :)
        real(dp), allocatable :: amplitudes(:)
        integer, allocatable :: set_flag(:)

        integer :: n_ops = 0

    contains

        procedure :: get_ex
        procedure :: get_ind
        procedure :: get_amp => get_amp_ind

    end type cc_amplitude

    ! and make a global cc_ops
    type(cc_amplitude), allocatable :: cc_ops(:), disc_cc_ops(:), all_cc_ops(:)

    logical :: t_store_disc_ops = .true.

    integer, allocatable :: ind_matrix_singles(:, :)
    integer, allocatable :: ind_matrix_doubles(:, :)

    integer :: n_virt_pairs = 0, nvirt = 0

    integer, allocatable :: elec_ind_mat(:, :), orb_ind_mat(:, :)

contains

    subroutine print_cc_amplitudes(hash_table, hash_size)
        ! routine to print out the cc-amplitudes to check there
        ! magnitude
        use util_mod, only: get_free_unit, get_unique_filename
        type(cc_hash), pointer, intent(in), optional :: hash_table(:)
        integer, intent(in), optional :: hash_size
        character(*), parameter :: this_routine = "print_cc_amplitudes"
        integer :: i, j, iunit
        character(12) :: filename
        character(1) :: x1
        type(cc_hash), pointer :: temp_node

        if (.not. allocated(cc_ops)) then
            print *, "cc amplitudes not yet allocated! cant print!"
            return
        end if

        if (present(hash_table)) then
            ASSERT(present(hash_size))

            iunit = get_free_unit()
            call get_unique_filename('cc_amps_4_wf', .true., .true., 1, &
                                     filename)

            open(iunit, file=filename, status='unknown')

            ! i have to do the quads special since it is stored in a hash..
            print *, "hash size_wf: ", hash_size
            do i = 1, hash_size
                temp_node => hash_table(i)

                if (temp_node%found) then
                    write(iunit, '(f15.7)') abs(temp_node%amp)

                end if

                do while (associated(temp_node%next))
                    if (temp_node%next%found) then
                        write(iunit, '(f15.7)') abs(temp_node%next%amp)
                    end if

                    temp_node => temp_node%next

                end do
            end do
            close(iunit)

            return
        end if

        if (iProcIndex == root) then
            do i = 1, 3
                ! open 4 files to print the cc-amps
                iunit = get_free_unit()
                write(x1, '(I1)') i
                call get_unique_filename('cc_amps_'//trim(x1), .true., .true., 1, &
                                         filename)
                open(iunit, file=filename, status='unknown')

                do j = 1, cc_ops(i)%n_ops
                    if (abs(cc_ops(i)%get_amp(j)) < EPS) cycle
                    write(iunit, '(f15.7)') abs(cc_ops(i)%get_amp(j))
                end do

                close(iunit)

            end do

            iunit = get_free_unit()
            call get_unique_filename('cc_amps_4', .true., .true., 1, &
                                     filename)

            open(iunit, file=filename, status='unknown')

            ! i have to do the quads special since it is stored in a hash..
            print *, "hash size: ", quad_hash_size
            do i = 1, quad_hash_size
                temp_node => quad_hash(i)

                if (temp_node%found) then
                    write(iunit, '(f15.7)') abs(temp_node%amp)

                end if

                do while (associated(temp_node%next))
                    if (temp_node%next%found) then
                        write(iunit, '(f15.7)') abs(temp_node%next%amp)
                    end if

                    temp_node => temp_node%next

                end do
            end do
            close(iunit)
        end if

    end subroutine print_cc_amplitudes

    subroutine calc_cc_quad_norm(hash_table, hash_size)
        ! need specific routine to calculate the quads norm, since it is
        ! stored in a hash-table format
        type(cc_hash), pointer, intent(in) :: hash_table(:)
        integer, intent(in) :: hash_size
        integer :: i, n_test
        type(cc_hash), pointer :: temp_node

        n_test = 0
        cc_amp_norm(:, 4) = 0.0_dp

        do i = 1, hash_size
            temp_node => hash_table(i)

            if (temp_node%found) then
                ! for now check if the loop over the hash table works:
                n_test = n_test + 1
                cc_amp_norm(1, 4) = cc_amp_norm(1, 4) + abs(temp_node%amp)
                cc_amp_norm(2, 4) = cc_amp_norm(2, 4) + temp_node%amp**2

            end if

            do while (associated(temp_node%next))
                temp_node => temp_node%next

                if (temp_node%found) then
                    n_test = n_test + 1

                    cc_amp_norm(1, 4) = cc_amp_norm(1, 4) + abs(temp_node%amp)
                    cc_amp_norm(2, 4) = cc_amp_norm(2, 4) + temp_node%amp**2

                end if
            end do
        end do

        root_print "checking if loop over hash table works as intented: "
        print *, "counted quads: ", n_test, "on Proc: ", iProcIndex
        print *, "L0 norm quads: ", cc_amp_norm(0, 4), "on proc: ", iProcIndex

        ! and apply square root to L2 Norm
        cc_amp_norm(2, 4) = sqrt(cc_amp_norm(2, 4))

    end subroutine calc_cc_quad_norm

    function cc_singles_factor() result(factor)
        ! this function should provide the correct factor to the
        ! cepa-shift for the singles.. if the correct variable are not
        ! yet set or sampled as 0 it should default to 1
        real(dp) :: factor

        real(dp) :: fac_triples, fac_doubles, weight

        weight = 1.0_dp / 4.0_dp

        ! the singles should be influenced by the triples and doubles..
        ! but this i have not figured out correctly..
        ! so for now return 1 always
!         factor = 1.0_dp
!         return

        if (cc_amp_norm(norm_comp, 3) < EPS) then
            fac_triples = 0.0_dp

            ! fix the weight in this case
            weight = 1.0_dp

        else
            ! for now just deal with the L^0 norm of the triples
            fac_triples = min(AllEXLEVEL_WNorm(norm_comp, 3, 1) &
                              / cc_amp_norm(norm_comp, 3), 1.0_dp)

        end if

        ! with the doubles i am scared that the estimated number of
        ! doubles could actually be lower then the sampled ones..
        if (cc_amp_norm(norm_comp, 2) < EPS) then

            fac_doubles = 0.0_dp

            weight = 0.0_dp

        else
            ! but 1 should be the maximum..
            fac_doubles = min(AllEXLEVEL_WNorm(0, 2, 1) &
                              / cc_amp_norm(norm_comp, 2), 1.0_dp)
        end if

        ! and then we have to combine the two factors with some weighting
        factor = 1.0_dp - (weight * fac_doubles + (1.0_dp - weight) * fac_triples)

    end function cc_singles_factor

    function cc_doubles_factor() result(factor)
        real(dp) :: factor

        real(dp) :: fac_triples, fac_quads, weight

        weight = 1.0_dp / 4.0_dp

        ! essentiall we would need the triples influence too..
        ! but Manu said this is not necessary.. hm.. lets try
        if (cc_amp_norm(norm_comp, 3) < EPS) then
            fac_triples = 0.0_dp
            weight = 0.0_dp
        else
            fac_triples = min(AllEXLEVEL_WNorm(0, 3, 1) / &
                              cc_amp_norm(norm_comp, 3), 1.0_dp)
        end if

        if (cc_amp_norm(norm_comp, 4) < EPS) then
            fac_quads = 0.0_dp
            weight = 1.0_dp
        else
            fac_quads = min(AllEXLEVEL_WNorm(0, 4, 1) / &
                            cc_amp_norm(norm_comp, 4), 1.0_dp)

        end if

        factor = 1.0_dp - (weight * fac_triples + (1.0_dp - weight) * fac_quads)

    end function cc_doubles_factor

    function cc_triples_factor() result(factor)
        real(dp) :: factor

        ! for now just set that to 1
        factor = 1.0_dp

        ! essentially we would need to have the influence of the
        ! quadrupels and the quintuples.. but those numbers are hard to
        ! estimate..
        ! todo

    end function cc_triples_factor

    function cc_quads_factor() result(factor)
        real(dp) :: factor

        ! see above
        factor = 1.0_dp

    end function cc_quads_factor

    subroutine setup_ind_matrix_doubles()
        character(*), parameter :: this_routine = "setup_ind_matrix_doubles"

        integer :: n_elec_pairs, i, j, a, b, elec_i, elec_j
        integer :: ij, ab, orb_a, orb_b, k
        logical :: t_par
        integer(n_int) :: temp_ilut(0:nifd)

        ASSERT(allocated(projedet))
        ASSERT(allocated(iLutRef))

        nvirt = nbasis - nel
        n_virt_pairs = nvirt * (nvirt - 1) / 2
        n_elec_pairs = nel * (nel - 1) / 2

        call setup_virtual_mask()

        call setup_elec_ind_mat()
        call setup_orb_ind_mat()

        ! encode only the possible excitations from the closed shell
        ! reference!
        ! i < j; a < b and (i,j) /= (a,b)

        ASSERT(n_elec_pairs == ElecPairs)

        allocate(ind_matrix_doubles(n_elec_pairs, n_virt_pairs))
        ind_matrix_doubles = 0

        temp_ilut = iLutRef(0:nifd, 1)
        k = 1

        ! i could also just loop over the reference determinant..
        do i = 1, nel
            elec_i = projedet(i, 1)
            do j = i + 1, nel
                elec_j = projedet(j, 1)

                t_par = same_spin(elec_i, elec_j)

                ! i have to convert (i,j) -> to a linear index
                ! and i definetly have to use the actual spin orbitals,
                ! since this is the quantity we have access to in the rest of
                ! the calculation
                ij = linear_elec_ind(elec_i, elec_j)
                ! maybe here i could to a check if ij = 0 ?

                ASSERT(ij > 0)
                ASSERT(ij <= ElecPairs)

                ! now i have to check if the orbitals fit and if the spin
                ! fits!
                ! use the virtual mask from back-spawn!
                if (t_par) then
                    do a = 1, Nvirt
                        orb_a = mask_virt_ni(a, 1)
                        do b = a + 1, Nvirt
                            orb_b = mask_virt_ni(b, 1)

                            if (same_spin(orb_a, orb_b) .and. same_spin(orb_a, elec_i)) then
                                ab = linear_orb_ind(orb_a, orb_b)
                                ASSERT(ab > 0)
                                ASSERT(ab <= n_virt_pairs)

                                ind_matrix_doubles(ij, ab) = k
                                k = k + 1

                            end if
                        end do
                    end do
                else
                    do a = 1, nvirt
                        orb_a = mask_virt_ni(a, 1)
                        do b = a + 1, nvirt
                            orb_b = mask_virt_ni(b, 1)
                            if (.not. same_spin(orb_a, orb_b)) then
                                ab = linear_orb_ind(orb_a, orb_b)

                                ASSERT(ab > 0)
                                ASSERT(ab <= n_virt_pairs)

                                ind_matrix_doubles(ij, ab) = k
                                k = k + 1
                            end if
                        end do
                    end do
                end if
            end do
        end do

    end subroutine setup_ind_matrix_doubles

    function linear_elec_ind(i, j) result(ind)
        integer, intent(in) :: i, j
        integer :: ind

        ! for a general orbital indexing i want to get a linear index
        ! for the electron pairs in the closed shell reference
        ! i could set up a matrix (nbasis, nbasis) for these pairs
        ! or i could make a search in the reference det at which position
        ! (i) and (j) are and use this then to encode the infomation..
        ! because i am lazy i think i will setup the matrix!
        ind = elec_ind_mat(i, j)

    end function linear_elec_ind

    function linear_orb_ind(a, b) result(ind)
        integer, intent(in) :: a, b
        integer :: ind

        ! see above!
        ind = orb_ind_mat(a, b)

    end function linear_orb_ind

    subroutine setup_elec_ind_mat
        character(*), parameter :: this_routine = "setup_elec_ind_mat"

        integer :: i, j, k, elec_i, elec_j
        if (allocated(elec_ind_mat)) deallocate(elec_ind_mat)
        ASSERT(allocated(projedet))

        allocate(elec_ind_mat(nbasis, nbasis))
        elec_ind_mat = 0

        k = 1

        do i = 1, nel
            elec_i = projedet(i, 1)
            do j = i + 1, nel
                elec_j = projedet(j, 1)

                elec_ind_mat(elec_i, elec_j) = k

                k = k + 1
            end do
        end do

        ASSERT(k - 1 == ElecPairs)

    end subroutine setup_elec_ind_mat

    subroutine setup_orb_ind_mat
        character(*), parameter :: this_routine = "setup_orb_ind_mat"

        integer :: i, j, k, orb_i, orb_j

        if (allocated(orb_ind_mat)) deallocate(orb_ind_mat)

        ASSERT(allocated(mask_virt_ni))

        allocate(orb_ind_mat(nbasis, nbasis))
        orb_ind_mat = 0

        k = 1
        do i = 1, Nvirt
            orb_i = mask_virt_ni(i, 1)
            do j = i + 1, nvirt
                orb_j = mask_virt_ni(j, 1)

                orb_ind_mat(orb_i, orb_j) = k
                k = k + 1
            end do
        end do

        ASSERT(k - 1 == n_virt_pairs)

    end subroutine setup_orb_ind_mat

    function get_amp_ind(this, ind) result(amp)
        ! write an amplitude getter, which gives 0 if the index does not
        ! fit
        class(cc_amplitude) :: this
        integer, intent(in) :: ind
        real(dp) :: amp

        if (ind == 0) then
            amp = 0.0_dp
        else
            amp = this%amplitudes(ind)
        end if

    end function get_amp_ind

    subroutine setup_ind_matrix_singles()
        character(*), parameter :: this_routine = "setup_ind_matrix_singles"

        integer :: i, j, k
        integer(n_int) :: temp_ilut(0:nifd)

        ASSERT(allocated(projedet))
        ASSERT(allocated(iLutRef))

        if (allocated(ind_matrix_singles)) deallocate(ind_matrix_singles)

        allocate(ind_matrix_singles(nbasis, nbasis))
        ind_matrix_singles = 0

        temp_ilut = iLutRef(0:nifd, 1)
        k = 1
        do i = 1, nbasis
            if (IsOcc(temp_ilut, i)) then
                ! the i is an electron in the reference!
                ! and for each electron index the virtuals
                do j = 1, nbasis
                    ! and i should and can specify the spin here!
                    if (IsNotOcc(temp_ilut, j) .and. same_spin(i, j)) then
                        ind_matrix_singles(i, j) = k
                        k = k + 1
                    end if
                end do
            end if
        end do

    end subroutine setup_ind_matrix_singles

    ! and now go to the routines to calculate the number of triples and
    ! quadrupels
    subroutine init_cc_amplitudes

        integer :: n_excits(cc_order)

        root_print "------ test on the cc ------- "

        n_excits = calc_number_of_excitations(nOccAlpha, nOccBeta, cc_order, &
                                              nbasis / 2)

        root_print "total number of possible excitations: "
        root_print n_excits

        call setup_ind_matrix_singles()
        call setup_ind_matrix_doubles()

        ! should i do this on each processors independently and then gather
        ! or should i gather first?
        ! if i want to calculate the T_2^2 i would need all of the
        ! doubles i guess.. so i can fill them independently, but should
        ! then communicat..
        call fill_cc_amplitudes(n_excits)
        ! to calculate all the possible cc-amps i need to communicate the
        ! operators.
        call communicate_cc_amps(n_excits)
        call calc_n_triples()
        if (t_store_hash_quadrupels) then
            quad_hash_size = n_doubles * (n_doubles - 1) / 2
            call init_cc_hash(quad_hash, quad_hash_size)
        end if
        call calc_n_quads()
        call calc_cc_quad_norm(quad_hash, quad_hash_size)

        print *, "number of hash clashes: ", n_clashes, "on proc: ", iProcIndex

        ! now i have to gather the information..
        call MPISumAll(n_singles, all_n_singles)
        call MPISumAll(n_doubles, all_n_doubles)
        call MPISumAll(est_triples, all_est_triples)
        ! can it do matrices:
        call MPISumAll(cc_amp_norm, all_cc_amp_norm)

        if (iProcIndex == root) then
            print *, " ---------- Singles -----------------"
            print *, "sampled singles in FCIQMC wavefunction: ", all_n_singles
            print *, "L0 Norm of singles in FCIQMC wavefunction: ", AllEXLEVEL_WNorm(0, 1, 1)
            print *, "L1 Norm of singles in FCIQMC wavefunction: ", AllEXLEVEL_WNorm(1, 1, 1)
            print *, "L2 Norm of singles in FCIQMC wavefunction: ", AllEXLEVEL_WNorm(2, 1, 1)

            print *, "L0 Norm of singles in CC-wavefunction: ", all_cc_amp_norm(0, 1)
            print *, "L1 Norm of singles in CC-wavefunction: ", all_cc_amp_norm(1, 1)
            print *, "L2 Norm of singles in CC-wavefunction: ", all_cc_amp_norm(2, 1)

            print *, ""
            print *, " ------------ Doubles -----------------"
            print *, "sampled doubles in FCIQMC wavefunction: ", all_n_doubles
            print *, "L0 Norm of doubles in FCIQMC wavefunction: ", AllEXLEVEL_WNorm(0, 2, 1)
            print *, "L1 Norm of doubles in FCIQMC wavefunction: ", AllEXLEVEL_WNorm(1, 2, 1)
            print *, "L2 Norm of doubles in FCIQMC wavefunction: ", AllEXLEVEL_WNorm(2, 2, 1)

            print *, "L0 Norm of doubles in CC-wavefunction: ", all_cc_amp_norm(0, 2)
            print *, "L1 Norm of doubles in CC-wavefunction: ", all_cc_amp_norm(1, 2)
            print *, "L2 Norm of doubles in CC-wavefunction: ", all_cc_amp_norm(2, 2)

            print *, ""
            print *, " ------------ Triples ------------------ "
            print *, "est. triples from T1^2 and T1*T2 ", all_est_triples
            print *, "L0 Norm of triples in FCIQMC wavefunction: ", AllEXLEVEL_WNorm(0, 3, 1)
            print *, "L1 Norm of triples in FCIQMC wavefunction: ", AllEXLEVEL_WNorm(1, 3, 1)
            print *, "L2 Norm of triples in FCIQMC wavefunction: ", AllEXLEVEL_WNorm(2, 3, 1)

            print *, "L0 Norm of triples in CC-wavefunction: ", all_cc_amp_norm(0, 3)
            print *, "L1 Norm of triples in CC-wavefunction: ", all_cc_amp_norm(1, 3)
            print *, "L2 Norm of triples in CC-wavefunction: ", all_cc_amp_norm(2, 3)

            print *, ""
            print *, " ----------- Quadruples ----------------- "
            print *, "L0 Norm of quadrupels in FCIQMC wavefunction: ", AllEXLEVEL_WNorm(0, 4, 1)
            print *, "L1 Norm of quadrupels in FCIQMC wavefunction: ", AllEXLEVEL_WNorm(1, 4, 1)
            print *, "L2 Norm of quadrupels in FCIQMC wavefunction: ", AllEXLEVEL_WNorm(2, 4, 1)

            print *, "estimated quadrupels from T2^2: "
            print *, "L0 Norm of quadrupels in CC-wavefunction: ", &
                all_cc_amp_norm(0, 4)
            print *, "L1 Norm of quadrupels in CC-wavefunction: ", &
                all_cc_amp_norm(1, 4)
            print *, "L2 Norm of quadrupels in CC-wavefunction: ", &
                all_cc_amp_norm(2, 4)

            print *, ""
        end if

    end subroutine init_cc_amplitudes

    subroutine communicate_cc_amps(n_excits)
        integer, intent(in) :: n_excits(cc_order)

        integer :: i, j
        ! only deal with systems with no single excitations, so we just
        ! have to communicate the doubles!
        allocate(all_cc_ops(2))

        do i = 1, 2
            all_cc_ops(i)%order = i
            allocate(all_cc_ops(i)%amplitudes(n_excits(i)))
            allocate(all_cc_ops(i)%operators(n_excits(i), 2, i))
            allocate(all_cc_ops(i)%set_flag(n_excits(i)))

            all_cc_ops(i)%operators = 0
            all_cc_ops(i)%amplitudes = 0.0_dp
            all_cc_ops(i)%set_flag = 0

            all_cc_ops(i)%n_ops = n_excits(i)

            call MPIReduce(cc_ops(i)%amplitudes, MPI_SUM, all_cc_ops(i)%amplitudes)
            do j = 1, i
                call MPIReduce(cc_ops(i)%operators(:, :, j), MPI_SUM, all_cc_ops(i)%operators(:, :, j))
            end do

            call MPIReduce(cc_ops(i)%set_flag, MPI_SUM, all_cc_ops(i)%set_flag)
        end do

        ! i just want to store up to triples or? but for these special
        ! case only consider systems with no single excitations..

    end subroutine communicate_cc_amps

    subroutine init_cc_hash(hash_table, hash_table_size)
        ! routine to setup the hash for the quadrupels
        type(cc_hash), pointer, intent(inout) :: hash_table(:)
        integer, intent(in) :: hash_table_size

        integer :: i
        ! estimate the hash table size with num_doubles^2 for now..
        ! although this number can be reduced substantially i guess

        allocate(hash_table(hash_table_size))

        do i = 1, hash_table_size
            hash_table(i)%found = .false.
            hash_table(i)%amp = 0.0_dp
            nullify (hash_table(i)%next)
        end do

    end subroutine init_cc_hash

    subroutine calc_n_triples()
        ! this routine calculates the number of "important" triples
        ! from the samples singles and doubles amplitudes:
        ! essentiall calulating T1 * T2

        integer :: i, j, ia(2, 1), jk_cd(2, 2)

        do i = 1, cc_ops(1)%n_ops
            ! for each t_i^a i have to check if a double excitation is
            ! possible on top of it..
            if (.not. cc_ops(1)%set_flag(i) == 1) cycle

            ia = cc_ops(1)%get_ex(i)

            do j = 1, cc_ops(2)%n_ops

                ! i also have to check if the amplitute was set
                if (.not. cc_ops(2)%set_flag(j) == 1) cycle

                jk_cd = cc_ops(2)%get_ex(j)

                ! check if the operators fit..
                if (any(ia(1, 1) == jk_cd(1, :)) .or. any(ia(2, 1) == jk_cd(2, :))) then
                    cycle
                end if
                ! if it fits increase the triples counter:
                est_triples = est_triples + 1

            end do
        end do

    end subroutine calc_n_triples

    subroutine calc_n_quads
        ! this routines calculates the number of "important" quadrupels
        ! calculating T2*T2 essentially
        ! todo: do i double count here? and how do i solve this issue?
        ! since different excitation generators can lead to the same
        ! quadruple excitation or?
        ! yes there are worst case 18 different permutations leading to the
        ! same quadruple excitation!
        ! so either divide by this number!
        ! or, for not too big numbers of double excitations create a
        ! specific hash table for the quadruple indexing out of the
        ! 8 indices i,j,k,l and a,b,c,d ..
        ! which would also allow me to store the estimated quadrupels
        ! magnitute(when ignoring the influence of the T1 and T3 operators!)
        ! but when i store the tripels also i could also estimate the
        ! quadrupels "exactly" .. which would be quite nice actually!
        ! todo: hash-table or dividing by 18..

        integer :: i, j, ij_ab(2, 2), kl_cd(2, 2), ijab_klcd(8), hash_ind
        integer(n_int) :: temp_int(0:nifd)
        real(dp) :: phase, amp
        logical :: t_found
        integer :: ijk_abc(2, 3), l_d(2, 1)

        print *, "test doubles: "
        do i = 1, cc_ops(2)%n_ops
            if (.not. cc_ops(2)%set_flag(i) == 1) cycle
            ij_ab = cc_ops(2)%get_ex(i)
            print *, "i:", i, "|", ij_ab(1, :), "->", ij_ab(2, :)
        end do

        if (t_store_hash_quadrupels) then
            do i = 1, cc_ops(2)%n_ops
                if (cc_ops(2)%set_flag(i) == 1) then

                    ij_ab = cc_ops(2)%get_ex(i)

                    do j = i + 1, cc_ops(2)%n_ops

                        if (cc_ops(2)%set_flag(j) == 1) then

                            kl_cd = cc_ops(2)%get_ex(j)

                            ! can i just: any(a==b)
                            if (unique_quad_ind_2_2(ij_ab, kl_cd)) then

                                ! i also want to store the amplitudes
                                ! i gues..
                                ! here i have to first order the indices
                                ! so i < j < k < l and a < b < c < d
                                ! and find the phase factor between them
                                call order_quad_indices_2_2(ij_ab, kl_cd, phase, &
                                                            ijab_klcd)

                                amp = phase * cc_ops(2)%get_amp(i) * cc_ops(2)%get_amp(j)

                                ! and i need a unique quantity associated with
                                ! these 8 orbital -> just set all the
                                ! corresponding orbitals
                                temp_int = apply_excit_ops(ilutRef(:, 1), &
                                                           [ij_ab(1, :), kl_cd(1, :)], [ij_ab(2, :), kl_cd(2, :)])

                                call cc_hash_look_up(ijab_klcd, temp_int, quad_hash, &
                                                     hash_ind, t_found)

                                ! it it is found in the hash-table i want to
                                ! update the amplitude
                                if (t_found) then
                                    print *, "repeated doubles: "
                                    print *, "i: ", i, "|", ij_ab
                                    print *, "j: ", j, "|", kl_cd

                                    call cc_hash_update(quad_hash, hash_ind, &
                                                        temp_int, amp)
                                else
                                    ! otherwise i want to add the entry
                                    print *, "initial doubles: "
                                    print *, "i: ", i, "|", ij_ab
                                    print *, "j: ", j, "|", kl_cd
                                    call cc_hash_add(quad_hash, hash_ind, &
                                                     temp_int, amp)

                                    ! and only in this case we found a new
                                    ! quadruple coefficient
                                    cc_amp_norm(0, 4) = cc_amp_norm(0, 4) + 1
                                end if
                            end if
                        end if
                    end do
                end if
            end do

            print *, "cc-4 L0 norm after T2^2", cc_amp_norm(0, 4), "on proc: ", &
                iProcIndex

            ! do i also want to do the T1*T3, T1^2*T2 T1^4
            do i = 1, cc_ops(3)%n_ops
                ijk_abc = cc_ops(3)%get_ex(i)

                do j = 1, cc_ops(1)%n_ops
                    if (cc_ops(1)%set_flag(j) == 1) then
                        l_d = cc_ops(1)%get_ex(j)

                        if (unique_quad_ind_3_1(ijk_abc, l_d)) then

                            call order_quad_indices_3_1(ijk_abc, l_d, phase, ijab_klcd)

                            amp = phase * cc_ops(3)%get_amp(i) * cc_ops(1)%get_amp(j)

                            temp_int = apply_excit_ops(ilutRef, &
                                                       [ijk_abc(1, :), l_d(1, 1)], [ijk_abc(2, :), l_d(2, 1)])

                            call cc_hash_look_up(ijab_klcd, temp_int, quad_hash, &
                                                 hash_ind, t_found)

                            ! it it is found in the hash-table i want to
                            ! update the amplitude
                            if (t_found) then
                                call cc_hash_update(quad_hash, hash_ind, &
                                                    temp_int, amp)
                            else
                                ! otherwise i want to add the entry
                                call cc_hash_add(quad_hash, hash_ind, &
                                                 temp_int, amp)

                                ! and only in this case we found a new
                                ! quadruple coefficient
                                cc_amp_norm(0, 4) = cc_amp_norm(0, 4) + 1
                            end if

                        end if
                    end if
                end do
            end do
            print *, "cc-4 L0 norm after T3*T1", cc_amp_norm(0, 4), "on proc: ", &
                iProcIndex

        else
            do i = 1, cc_ops(2)%n_ops
                if (cc_ops(2)%set_flag(i) == 1) then

                    ij_ab = cc_ops(2)%get_ex(i)

                    do j = i + 1, cc_ops(2)%n_ops

                        if (cc_ops(2)%set_flag(j) == 1) then

                            kl_cd = cc_ops(2)%get_ex(j)

                            ! can i just: any(a==b)
                            if (unique_quad_ind_2_2(ij_ab, kl_cd)) then

                                cc_amp_norm(0, 4) = cc_amp_norm(0, 4) + 1

                            end if
                        end if
                    end do
                end if
            end do
        end if

    end subroutine calc_n_quads

    subroutine cc_hash_look_up(ind, tgt, hash_table, hash_val, t_found)
        integer, intent(in) :: ind(:)
        integer(n_int), intent(in) :: tgt(0:nifd)
        type(cc_hash), pointer :: hash_table(:)
        integer, intent(out) :: hash_val
        logical, intent(out) :: t_found

        type(cc_hash), pointer :: temp_node

        t_found = .false.

        hash_val = FindWalkerHash(ind, size(hash_table))

        temp_node => hash_table(hash_val)

        if (temp_node%found) then
            do while (associated(temp_node))
                if (all(tgt == temp_node%ind)) then
                    t_found = .true.
                    exit
                end if
                temp_node => temp_node%next
            end do
        end if

        nullify (temp_node)

    end subroutine cc_hash_look_up

    subroutine cc_hash_update(hash_table, hash_val, tgt, amp)
        ! routine to update the amplitude of a found entry
        type(cc_hash), pointer, intent(inout) :: hash_table(:)
        integer, intent(in) :: hash_val
        integer(n_int), intent(in) :: tgt(:)
        real(dp), intent(in) :: amp
        character(*), parameter :: this_routine = "cc_hash_update"

        type(cc_hash), pointer :: temp_node
        logical :: found

        found = .false.

        temp_node => hash_table(hash_val)

        do while (associated(temp_node))
            if (all(temp_node%ind == tgt)) then
                ! this is the correct entry!
                found = .true.
                temp_node%amp = temp_node%amp + amp
                exit
            end if
            temp_node => temp_node%next
        end do

        ASSERT(found)

    end subroutine cc_hash_update

    subroutine cc_hash_add(hash_table, hash_val, tgt, amp)
        ! this is a routine to add a hash table entry
        ! it means the entry was not there yet or has to be added to
        ! the linked list
        type(cc_hash), pointer, intent(inout) :: hash_table(:)
        integer, intent(in) :: hash_val
        integer(n_int), intent(in) :: tgt(:)
        real(dp), intent(in) :: amp

        type(cc_hash), pointer :: temp_node

        temp_node => hash_table(hash_val)

        if (.not. temp_node%found) then
            ! this means this entry is empty so we cann fill it here
            temp_node%found = .true.
            allocate(temp_node%ind(0:nifd))
            temp_node%ind = tgt
            temp_node%amp = amp
        else
            ! loop through the linked list
            do while (associated(temp_node%next))
                temp_node => temp_node%next
            end do
            allocate(temp_node%next)
            nullify (temp_node%next%next)
            temp_node%next%found = .true.
            allocate(temp_node%next%ind(0:nifd))
            temp_node%next%ind = tgt
            temp_node%next%amp = amp

            ! for testing cound the number of clashes:
            n_clashes = n_clashes + 1
        end if

        nullify (temp_node)

    end subroutine cc_hash_add

    subroutine order_quad_indices_3_1(ijk_abc, l_d, phase, ijab_klcd)
        integer, intent(inout) :: ijk_abc(2, 3), l_d(2, 1)
        real(dp), intent(out) :: phase
        integer, intent(out) :: ijab_klcd(8)
        character(*), parameter :: this_routine = "order_quad_indices_3_1"

        call stop_all(this_routine, "TODO")
        unused_var(ijk_abc)
        unused_var(l_d)
        phase = 0.0_dp
        ijab_klcd = 0

    end subroutine order_quad_indices_3_1

    subroutine order_quad_indices_2_2(ij_ab, kl_cd, phase, ijab_klcd)
        integer, intent(inout) :: ij_ab(2, 2), kl_cd(2, 2)
        real(dp), intent(out) :: phase
        integer, intent(out) :: ijab_klcd(8)
        character(*), parameter :: this_routine = "order_quad_indices_2_2"

        integer :: ij(2), ab(2), kl(2), cd(2)
        integer :: n
        ! the correctly order t_ij^ab * t_kl^cd has the sign convention -1!
        ! and we can be sure that the inputed (ij),(ab) etc. are ordered i < j
        ! within the T2 operators, since only those get stored by convention
        ! and also that i < k and a < c

        ij = ij_ab(1, :)
        ab = ij_ab(2, :)
        kl = kl_cd(1, :)
        cd = kl_cd(2, :)

        n = 1

        ! assert i < k and a < c
        ASSERT(ij_ab(1, 1) < kl_cd(1, 1))
        ASSERT(ij(1) < ij(2))
        ASSERT(cd(1) < cd(2))
        ASSERT(ab(1) < ab(2))
        ASSERT(cd(1) < cd(2))

        ! with this: i am not so sure: a < c ..
        ! ok i can not be sure that a < c! so i have to sort more here!

        ! sort the electrons:
        ! if j < k everything is fine and in order!
        if (ij_ab(1, 2) < kl_cd(1, 1)) then
            ! everything is fine..

        else if (ij_ab(1, 2) < kl_cd(1, 2)) then
            ! if j < l then i have to switch j and k
            call swap(ij(2), kl(1))
            n = n + 1

        else
            ! this means j > l so we have to swap j -> l and then k -> l
            call swap(ij(2), kl(2))
            call swap(ij(2), kl(1))
            ! wich leaves the phase unchanged

        end if

        ! and check if we did everything correctly:
        ASSERT(ij(1) < ij(2))
        ASSERT(kl(1) < kl(2))
        ASSERT(ij(2) < kl(1))

        ! i am not so sure any more if a < c anymore
        ! now sort the orbitals but here a > c is also possible!
        if (ij_ab(2, 1) < kl_cd(2, 1)) then
            ! so a < c
            if (ij_ab(2, 2) < kl_cd(2, 1)) then
                ! b < c so in this case nothing has to be done!

            else if (ij_ab(2, 2) < kl_cd(2, 2)) then
                ! b < d: so b and c have to be swapped!
                call swap(ab(2), cd(1))
                n = n + 1

            else
                ! b > d: so switch b -> d and d -> c
                call swap(ab(2), cd(2))
                call swap(ab(2), cd(1))
                ! and this does not change the phase
            end if
        else
            ! so this means a > c which implies b > c also!
            ! so if b < d we know everything
            if (ij_ab(2, 2) < kl_cd(2, 2)) then
                ! then c < a < b < d
                call swap(ab(1), cd(1))
                call swap(ab(2), cd(1))
                ! so no phase factor!

            else if (ij_ab(2, 1) > kl_cd(2, 2)) then
                ! here b and a > d
                ! so c < d < a < b
                call swap(ab(1), cd(1))
                call swap(ab(2), cd(2))
                ! also here no phase factor!

            else
                ! this means c < a < d < b
                call swap(ab(1), ab(2))
                ! ba, cd
                call swap(cd(1), cd(2))
                ! ba, dc
                call swap(ab(1), cd(2))
                ! cadb
                ! here we get a phase
                n = n + 1
            end if
        end if

        ASSERT(ab(1) < ab(2))
        ASSERT(cd(1) < cd(2))
        ASSERT(ab(2) < cd(1))

        ! now switch the indices
        ij_ab(1, :) = ij
        ij_ab(2, :) = ab
        kl_cd(1, :) = kl
        kl_cd(2, :) = cd

        ! and also store a linear index
        ijab_klcd = [ij, ab, kl, cd]

        ! and calculate the appropriate phase
        phase = (-1.0_dp)**n

    end subroutine order_quad_indices_2_2

    logical function unique_quad_ind_3_1(ijk_abc, l_d)
        integer, intent(in) :: ijk_abc(2, 3), l_d(2, 1)

        unique_quad_ind_3_1 = all(l_d(1, 1) /= ijk_abc(1, :) .and. &
                                  l_d(2, 1) /= ijk_abc(2, :))

    end function unique_quad_ind_3_1

    logical function unique_quad_ind_2_2(ij_ab, kl_cd)
        integer, intent(in) :: ij_ab(2, 2), kl_cd(2, 2)

        unique_quad_ind_2_2 = all(ij_ab(1, 1) /= kl_cd(1, :) .and. ij_ab(1, 2) /= kl_cd(1, :) &
                                  .and. ij_ab(2, 1) /= kl_cd(2, :) .and. ij_ab(2, 2) /= kl_cd(2, :))

    end function unique_quad_ind_2_2

    function apply_excit_ops(ilut_in, elecs, orbs) result(ilut_out)
        integer(n_int), intent(in) :: ilut_in(0:nifd)
        integer, intent(in) :: elecs(:), orbs(:)
        integer(n_int) :: ilut_out(0:nifd)

        integer :: i

        ilut_out = ilut_in

        do i = 1, size(elecs)
            associate(j => elecs(i), k => orbs(i))
                clr_orb(ilut_out, j)
                set_orb(ilut_out, k)
            end associate
        end do

    end function apply_excit_ops

    ! do it other way.. only store the possible non-zero cluster operators!
    subroutine fill_cc_amplitudes(n_excits)
        ! design decisions: since the singles are not so many in general
        ! and because i could need them to correct the doubles amplitudes
        ! store all of the possible ones! and encode them specifically through
        ! (i) and (a)
        integer, intent(in) :: n_excits(cc_order)

        character(*), parameter :: this_routine = "fill_cc_amplitudes"

        integer :: idet, ic, ex(2, cc_order), j, i
        integer :: ia, ib, ja, jb, ind

        integer :: a, b, elec_i, elec_j, orb_a, orb_b
        logical :: t_par, t_store_full_doubles = .true.

        real(dp) :: amp

        real(dp), allocatable :: temp_amps(:)
        integer, allocatable :: temp_ops(:, :, :)

        integer :: ab(2), ac(2), bc(2), ij(2), ik(2), jk(2)
        integer :: c, k, ij_ac, ij_bc, ij_ab, ik_bc, ik_ac, ik_ab
        integer :: jk_bc, jk_ac, jk_ab, elec_k, orb_c, jc, ka, kb, kc, n
        integer(n_int) :: temp_ilut(0:nifd)
        integer :: dummy_ind, dummy_hash, temp_nI(nel)
        logical :: tSuccess
        integer :: i_a(2, 1), j_b(2, 1)
        integer :: n_disc_1, n_disc_2
        integer :: quad_ind(8), hash_ind
        logical :: t_found
        integer(n_int) :: temp_int(0:nifd)

        real(dp) :: sign_tmp(lenof_sign)

        ! for this it is helpful to have an upper limit of the number of
        ! possible amplitudes, but just do it for the singles for now..

        allocate(cc_ops(cc_order))

        if (t_store_full_doubles .and. cc_order > 2 .and. t_store_hash_quadrupels) then
            ! also store the intermediate T1^2 and T2^2 in operator form to
            ! make my life easier
            allocate(disc_cc_ops(2))
            disc_cc_ops(1)%order = 2
            disc_cc_ops(2)%order = 4
        end if

        ! and do a nice initialization depending on the order
        do i = 1, cc_order
            cc_ops(i)%order = i
        end do

        allocate(cc_ops(1)%amplitudes(n_excits(1)))
        allocate(cc_ops(1)%operators(n_excits(1), 2, 1))
        allocate(cc_ops(1)%set_flag(n_excits(1)))

        cc_ops(1)%operators = 0
        cc_ops(1)%amplitudes = 0.0_dp
        cc_ops(1)%set_flag = 0

        cc_ops(1)%n_ops = n_excits(1)

        ! first figure out the number of double and fill the singles
        ! amplitudes!

        n_doubles = 0
        do idet = 1, int(totwalkers)
            ! for now also count the exact number of triples and quadrupels
            ic = FindBitExcitLevel(ilutRef(:, 1), CurrentDets(:, idet))
            call extract_sign(CurrentDets(:, idet), sign_tmp)

            ! for some strange reason there are zero elements stored..
            if (abs(sign_tmp(1)) > EPS) then
                select case (ic)
                case (1)
                    n_singles = n_singles + 1
                    call get_bit_excitmat(iLutRef(:, 1), CurrentDets(:, idet), ex, ic)
                    ind = cc_ops(1)%get_ind(ex(1, 1), ex(2, 1))

                    ! for now only do it for single runs.. do i need the normalising?
                    amp = sign_tmp(1)

                    cc_ops(1)%amplitudes(ind) = amp
                    cc_ops(1)%operators(ind, 1, 1) = ex(1, 1)
                    cc_ops(1)%operators(ind, 2, 1) = ex(2, 1)
                    cc_ops(1)%set_flag(ind) = 1
                    cc_amp_norm(0, 1) = cc_amp_norm(0, 1) + 1
                    cc_amp_norm(1, 1) = cc_amp_norm(1, 1) + abs(amp)
                    cc_amp_norm(2, 1) = cc_amp_norm(2, 1) + amp**2

                case (2)
                    n_doubles = n_doubles + 1

                case (3)
                    n_triples = n_triples + 1

                case (4)
                    n_quads = n_quads + 1

                end select
            end if
        end do

        cc_amp_norm(2, 1) = sqrt(cc_amp_norm(2, 1))

        print *, "direct sampled doubles: ", n_doubles, "on proc: ", iProcIndex

        if (allocated(cc_ops(2)%operators)) deallocate(cc_ops(2)%operators)
        if (allocated(cc_ops(2)%amplitudes)) deallocate(cc_ops(2)%amplitudes)

        ! here i have a choice to only store the non-zero contributions
        ! or encode them in the full-list to better access them later on.
        ! especially when we want to calculate the coupled-cluster
        ! triples and quadrupels..

        if (t_store_full_doubles) then
            allocate(cc_ops(2)%operators(n_excits(2), 2, 2))
            allocate(cc_ops(2)%amplitudes(n_excits(2)))
            allocate(cc_ops(2)%set_flag(n_excits(2)))
            cc_ops(2)%n_ops = n_excits(2)

            ! then also store the intermediate
            allocate(disc_cc_ops(1)%operators(n_excits(2), 2, 2))
            allocate(disc_cc_ops(1)%amplitudes(n_excits(2)))
            allocate(disc_cc_ops(1)%set_flag(n_excits(2)))
            disc_cc_ops(1)%n_ops = n_excits(2)

            disc_cc_ops(1)%operators = 0
            disc_cc_ops(1)%amplitudes = 0
            disc_cc_ops(1)%set_flag = 0

        else
            allocate(cc_ops(2)%operators(n_doubles, 2, 2))
            allocate(cc_ops(2)%amplitudes(n_doubles))
            allocate(cc_ops(2)%set_flag(n_doubles))

            cc_ops(2)%n_ops = n_doubles

            cc_ops(2)%operators = 0
            cc_ops(2)%amplitudes = 0.0_dp
            cc_ops(2)%set_flag = 0
        end if

        if (t_store_disc_ops) then
            n_disc_1 = 0
            n_disc_2 = 0
            do i = 1, cc_ops(1)%n_ops
                if (cc_ops(1)%set_flag(i) == 1) then
                    do j = i + 1, cc_ops(1)%n_ops
                        if (cc_ops(1)%set_flag(j) == 1) then
                            i_a = cc_ops(1)%get_ex(i)
                            j_b = cc_ops(1)%get_ex(j)

                            ! ensure ordering, they should be ordered by
                            ! electrons but electron index can still be the
                            ! same ..
                            if (i_a(1, 1) > j_b(1, 1)) then
                                call stop_all(this_routine, &
                                              "i did smth wrong with the ordering")
                            end if

                            if (i_a(1, 1) /= j_b(1, 1)) then
                                amp = cc_ops(1)%get_amp(i) * cc_ops(1)%get_amp(j)

                                if (i_a(2, 1) < j_b(2, 1)) then
                                    ind = cc_ops(2)%get_ind([i_a(1, 1), j_b(1, 1)], &
                                                            [i_a(2, 1), j_b(2, 1)])
                                    amp = -amp

                                else if (i_a(2, 1) > j_b(2, 1)) then
                                    ind = cc_ops(2)%get_ind([i_a(1, 1), j_b(1, 1)], &
                                                            [j_b(2, 1), i_a(2, 1)])

                                else
                                    ! orbitals are the same..
                                end if

                                ! did i already find that:
                                if (disc_cc_ops(1)%set_flag(ind) == 1) then
                                    disc_cc_ops(1)%amplitudes(ind) = &
                                        disc_cc_ops(1)%amplitudes(ind) + amp

                                    ! if it cancelled: remove the entry:
                                    if (abs(disc_cc_ops(1)%amplitudes(ind)) < EPS) then
                                        disc_cc_ops(1)%set_flag(ind) = 0
                                        disc_cc_ops(1)%operators(ind, :, :) = 0
                                        disc_cc_ops(1)%amplitudes(ind) = 0.0_dp
                                        ! and lower ofc..
                                        n_disc_1 = n_disc_1 - 1
                                    end if
                                else
                                    ! otherwise set it..
                                    disc_cc_ops(1)%amplitudes(ind) = amp
                                    disc_cc_ops(1)%set_flag(ind) = 1
                                    disc_cc_ops(1)%operators(ind, 1, :) = [i_a(1, 1), j_b(1, 1)]
                                    disc_cc_ops(1)%operators(ind, 2, :) = &
                                        [min(i_a(2, 1), j_b(2, 1)), max(i_a(2, 1), j_b(2, 1))]
                                    n_disc_1 = n_disc_1 + 1
                                end if
                            end if
                        end if
                    end do
                end if
            end do
        end if

        print *, "number of disconnected T1^2: ", n_disc_1, "on proc: ", iProcIndex
        cc_ops(2)%operators = 0
        cc_ops(2)%amplitudes = 0.0_dp
        cc_ops(2)%set_flag = 0

        do idet = 1, int(totwalkers)

            ic = FindBitExcitLevel(ilutRef(:, 1), CurrentDets(:, idet))

            if (ic == 2) then
                call extract_sign(CurrentDets(:, idet), sign_tmp)
                ! i hope this routine is not buggy..
                call get_bit_excitmat(iLutRef(:, 1), CurrentDets(:, idet), ex, ic)

                ia = cc_ops(1)%get_ind(ex(1, 1), ex(2, 1))
                ib = cc_ops(1)%get_ind(ex(1, 1), ex(2, 2))
                ja = cc_ops(1)%get_ind(ex(1, 2), ex(2, 1))
                jb = cc_ops(1)%get_ind(ex(1, 2), ex(2, 2))

                ! check first if the amp gets 0 due to the singles
                amp = sign_tmp(1) + &
                      cc_ops(1)%get_amp(ja) * cc_ops(1)%get_amp(ib) - &
                      cc_ops(1)%get_amp(ia) * cc_ops(1)%get_amp(jb)

                if (abs(amp) > EPS) then
                    if (t_store_full_doubles) then
                        ind = cc_ops(2)%get_ind(ex(1, :), ex(2, :))
                    else
                        ind = j
                        j = j + 1
                    end if

                    cc_ops(2)%operators(ind, 1, :) = ex(1, 1:2)
                    cc_ops(2)%operators(ind, 2, :) = ex(2, 1:2)
                    cc_ops(2)%amplitudes(ind) = amp
                    cc_ops(2)%set_flag(ind) = 1

                    ! i think i can already sum up the norms here..
                    ! since we have all the necessary contribs
                    cc_amp_norm(0, 2) = cc_amp_norm(0, 2) + 1
                    cc_amp_norm(1, 2) = cc_amp_norm(1, 2) + abs(amp)
                    cc_amp_norm(2, 2) = cc_amp_norm(2, 2) + amp**2

                end if
            end if
        end do

        ! i just realise that i have to run over all the possible double
        ! excitations, since t_ij^ab = c_ij^ab - t_i^a... could be
        ! non-zero, just from the single contribution..
        if (t_store_full_doubles) then
            do i = 1, nel
                elec_i = projedet(i, 1)
                do j = i + 1, nel
                    elec_j = projedet(j, 1)
                    t_par = same_spin(elec_i, elec_j)

                    if (t_par) then
                        do a = 1, nvirt
                            orb_a = mask_virt_ni(a, 1)
                            do b = a + 1, nvirt
                                orb_b = mask_virt_ni(b, 1)

                                if (same_spin(orb_a, orb_b) .and. &
                                    same_spin(orb_a, elec_i)) then

                                    ind = cc_ops(2)%get_ind([elec_i, elec_j], &
                                                            [orb_a, orb_b])

                                    ASSERT(elec_i < elec_j)
                                    ASSERT(orb_a < orb_b)

                                    ! do the disc op storage here!
                                    ia = cc_ops(1)%get_ind([elec_i], [orb_a])
                                    ib = cc_ops(1)%get_ind([elec_i], [orb_b])
                                    ja = cc_ops(1)%get_ind([elec_j], [orb_a])
                                    jb = cc_ops(1)%get_ind([elec_j], [orb_b])

                                    amp = cc_ops(1)%get_amp(ja) * cc_ops(1)%get_amp(ib) - &
                                          cc_ops(1)%get_amp(ia) * cc_ops(1)%get_amp(jb)

                                    if (abs(amp) > EPS) then
                                        if (.not. cc_ops(2)%set_flag(ind) == 1) then
                                            print *, "para: ", elec_i, elec_j, orb_a, orb_b
                                            cc_ops(2)%operators(ind, 1, :) = [elec_i, elec_j]
                                            cc_ops(2)%operators(ind, 2, :) = [orb_a, orb_b]
                                            cc_ops(2)%amplitudes(ind) = amp
                                            cc_ops(2)%set_flag(ind) = 1

                                            n_doubles = n_doubles + 1

                                            cc_amp_norm(0, 2) = cc_amp_norm(0, 2) + 1
                                            cc_amp_norm(1, 2) = cc_amp_norm(1, 2) + abs(amp)
                                            cc_amp_norm(2, 2) = cc_amp_norm(2, 2) + amp**2
                                        end if
                                        if (t_store_disc_ops) then

                                            disc_cc_ops(1)%operators(ind, 1, :) = [elec_i, elec_j]
                                            disc_cc_ops(1)%operators(ind, 2, :) = [orb_a, orb_b]
                                            disc_cc_ops(1)%amplitudes(ind) = amp
                                            disc_cc_ops(1)%set_flag(ind) = 1

                                            n_disc_2 = n_disc_2 + 1

                                        end if
                                    end if
                                end if
                            end do
                        end do
                    else
                        do a = 1, nvirt
                            orb_a = mask_virt_ni(a, 1)
                            do b = a + 1, nvirt
                                orb_b = mask_virt_ni(b, 1)

                                if (.not. same_spin(orb_a, orb_b)) then

                                    ind = cc_ops(2)%get_ind([elec_i, elec_j], &
                                                            [orb_a, orb_b])

                                    ASSERT(elec_i < elec_j)
                                    ASSERT(orb_a < orb_b)

                                    ia = cc_ops(1)%get_ind([elec_i], [orb_a])
                                    ib = cc_ops(1)%get_ind([elec_i], [orb_b])
                                    ja = cc_ops(1)%get_ind([elec_j], [orb_a])
                                    jb = cc_ops(1)%get_ind([elec_j], [orb_b])

                                    amp = cc_ops(1)%get_amp(ja) * cc_ops(1)%get_amp(ib) - &
                                          cc_ops(1)%get_amp(ia) * cc_ops(1)%get_amp(jb)

                                    if (abs(amp) > EPS) then

                                        if (.not. cc_ops(2)%set_flag(ind) == 1) then

                                            print *, "anti: ", elec_i, elec_j, orb_a, orb_b
                                            cc_ops(2)%operators(ind, 1, :) = [elec_i, elec_j]
                                            cc_ops(2)%operators(ind, 2, :) = [orb_a, orb_b]
                                            cc_ops(2)%amplitudes(ind) = amp
                                            cc_ops(2)%set_flag(ind) = 1

                                            n_doubles = n_doubles + 1

                                            cc_amp_norm(0, 2) = cc_amp_norm(0, 2) + 1
                                            cc_amp_norm(1, 2) = cc_amp_norm(1, 2) + abs(amp)
                                            cc_amp_norm(2, 2) = cc_amp_norm(2, 2) + amp**2
                                        end if
                                        if (t_store_disc_ops) then

                                            disc_cc_ops(1)%operators(ind, 1, :) = [elec_i, elec_j]
                                            disc_cc_ops(1)%operators(ind, 2, :) = [orb_a, orb_b]
                                            disc_cc_ops(1)%amplitudes(ind) = amp
                                            disc_cc_ops(1)%set_flag(ind) = 1

                                            n_disc_2 = n_disc_2 + 1
                                        end if

                                    end if
                                end if
                            end do
                        end do
                    end if
                end do
            end do
!             end if
        end if

        print *, "disconnecte T1^2 calc. differently: ", n_disc_2, "on proc: ", &
            iProcIndex

        cc_amp_norm(2, 2) = sqrt(cc_amp_norm(2, 2))

        print *, "doubles after singles contribution: ", n_doubles, "on proc:", &
            iProcIndex

        ! and should i also do the triples.. just to check maybe?
        ! but here i definetly only want to store the non-zero
        ! contributions.. but.. i actually would need to run over
        ! all the possible ones also since the contribution could be
        ! approximated by T1*T2 or T1^3..

        ! i could loop.. but i do not want to encode all of them!
        ! but essentially i have to first allocate a vector of all possible
        ! ones to be sure..
        if (cc_order >= 3) then
            allocate(temp_amps(n_excits(3)))
            allocate(temp_ops(n_excits(3), 2, 3))

            n = 1

            ! first analyze the wavefunction:
            do idet = 1, int(totwalkers)
                ic = FindBitExcitLevel(iLutRef(:, 1), CurrentDets(:, idet))

                if (ic == 3) then

                    call extract_sign(CurrentDets(:, idet), sign_tmp)
                    call get_bit_excitmat(iLutRef(:, 1), CurrentDets(:, idet), ex, ic)

                    ia = cc_ops(1)%get_ind(ex(1, 1), ex(2, 1))
                    ib = cc_ops(1)%get_ind(ex(1, 1), ex(2, 2))
                    ic = cc_ops(1)%get_ind(ex(1, 1), ex(2, 3))

                    ja = cc_ops(1)%get_ind(ex(1, 2), ex(2, 1))
                    jb = cc_ops(1)%get_ind(ex(1, 2), ex(2, 2))
                    jc = cc_ops(1)%get_ind(ex(1, 2), ex(2, 3))

                    ka = cc_ops(1)%get_ind(ex(1, 3), ex(2, 1))
                    kb = cc_ops(1)%get_ind(ex(1, 3), ex(2, 2))
                    kc = cc_ops(1)%get_ind(ex(1, 3), ex(2, 3))

                    ij = [ex(1, 1), ex(1, 2)]
                    ik = [ex(1, 1), ex(1, 3)]
                    jk = [ex(1, 2), ex(1, 3)]

                    ab = [ex(2, 1), ex(2, 2)]
                    ac = [ex(2, 1), ex(2, 3)]
                    bc = [ex(2, 2), ex(2, 3)]

                    jk_bc = cc_ops(2)%get_ind(jk, bc)
                    jk_ac = cc_ops(2)%get_ind(jk, ac)
                    jk_ab = cc_ops(2)%get_ind(jk, ab)

                    ik_bc = cc_ops(2)%get_ind(ik, bc)
                    ik_ac = cc_ops(2)%get_ind(ik, ac)
                    ik_ab = cc_ops(2)%get_ind(ik, ab)

                    ij_bc = cc_ops(2)%get_ind(ij, bc)
                    ij_ac = cc_ops(2)%get_ind(ij, ac)
                    ij_ab = cc_ops(2)%get_ind(ij, ab)

                    amp = sign_tmp(1) &
                          - cc_ops(1)%get_amp(ia) * cc_ops(2)%get_amp(jk_bc) &
                          + cc_ops(1)%get_amp(ib) * cc_ops(2)%get_amp(jk_ac) &
                          - cc_ops(1)%get_amp(ic) * cc_ops(2)%get_amp(jk_ab) &
                          + cc_ops(1)%get_amp(ja) * cc_ops(2)%get_amp(ik_bc) &
                          - cc_ops(1)%get_amp(jb) * cc_ops(2)%get_amp(ik_ac) &
                          + cc_ops(1)%get_amp(jc) * cc_ops(2)%get_amp(ik_ab) &
                          - cc_ops(1)%get_amp(ka) * cc_ops(2)%get_amp(ij_bc) &
                          + cc_ops(1)%get_amp(kb) * cc_ops(2)%get_amp(ij_ac) &
                          - cc_ops(1)%get_amp(kc) * cc_ops(2)%get_amp(ij_ab) &
                          - cc_ops(1)%get_amp(ia) * cc_ops(1)%get_amp(jb) * cc_ops(1)%get_amp(kc) &
                          + cc_ops(1)%get_amp(ia) * cc_ops(1)%get_amp(jc) * cc_ops(1)%get_amp(kb) &
                          + cc_ops(1)%get_amp(ib) * cc_ops(1)%get_amp(ja) * cc_ops(1)%get_amp(kc) &
                          - cc_ops(1)%get_amp(ib) * cc_ops(1)%get_amp(jc) * cc_ops(1)%get_amp(ka) &
                          - cc_ops(1)%get_amp(ic) * cc_ops(1)%get_amp(ja) * cc_ops(1)%get_amp(kb) &
                          + cc_ops(1)%get_amp(ic) * cc_ops(1)%get_amp(jb) * cc_ops(1)%get_amp(ka)

                    if (abs(amp) > EPS) then

                        temp_amps(n) = amp
                        temp_ops(n, :, :) = ex

                        n = n + 1

                        ! also here i can add the norms up already..
                        cc_amp_norm(0, 3) = cc_amp_norm(0, 3) + 1
                        cc_amp_norm(1, 3) = cc_amp_norm(1, 3) + abs(amp)
                        cc_amp_norm(2, 3) = cc_amp_norm(2, 3) + amp**2

                    end if
                end if
            end do
            print *, "direct sampled triples: ", n - 1, "on Proc: ", iProcIndex

            do i = 1, nel
                elec_i = projedet(i, 1)
                do j = i + 1, nel
                    elec_j = projedet(j, 1)
                    ij = [elec_i, elec_j]
                    do k = j + 1, nel
                        elec_k = projedet(k, 1)

                        ik = [elec_i, elec_k]
                        jk = [elec_j, elec_k]

                        do a = 1, nvirt
                            orb_a = mask_virt_ni(a, 1)

                            ia = cc_ops(1)%get_ind([elec_i], [orb_a])
                            ja = cc_ops(1)%get_ind([elec_j], [orb_a])
                            ka = cc_ops(1)%get_ind([elec_k], [orb_a])

                            do b = a + 1, nvirt
                                orb_b = mask_virt_ni(b, 1)
                                ab = [orb_a, orb_b]
                                ib = cc_ops(1)%get_ind([elec_i], [orb_b])
                                jb = cc_ops(1)%get_ind([elec_j], [orb_b])
                                kb = cc_ops(1)%get_ind([elec_k], [orb_b])

                                ij_ab = cc_ops(2)%get_ind(ij, ab)
                                ik_ab = cc_ops(2)%get_ind(ik, ab)
                                jk_ab = cc_ops(2)%get_ind(jk, ab)

                                do c = b + 1, nvirt
                                    orb_c = mask_virt_ni(c, 1)
                                    ac = [orb_a, orb_c]
                                    bc = [orb_b, orb_c]

                                    ic = cc_ops(1)%get_ind([elec_i], [orb_c])
                                    jc = cc_ops(1)%get_ind([elec_j], [orb_c])
                                    kc = cc_ops(1)%get_ind([elec_k], [orb_c])

                                    ij_ac = cc_ops(2)%get_ind(ij, ac)
                                    ik_ac = cc_ops(2)%get_ind(ik, ac)
                                    jk_ac = cc_ops(2)%get_ind(jk, ac)

                                    ij_bc = cc_ops(2)%get_ind(ij, bc)
                                    ik_bc = cc_ops(2)%get_ind(ik, bc)
                                    jk_bc = cc_ops(2)%get_ind(jk, bc)

                                    ! and now i have to check if the
                                    ! triple excitation (ijk|abc) is in the
                                    ! wallker list! use the hash-table!
                                    ! todo
                                    temp_ilut = apply_excit_ops(iLutRef, &
                                                                [elec_i, elec_j, elec_k], [orb_a, orb_b, orb_c])

                                    call decode_bit_det(temp_nI, temp_ilut)

                                    call hash_table_lookup(temp_nI, temp_ilut, &
                                                           nifd, HashIndex, CurrentDets, dummy_ind, &
                                                           dummy_hash, tSuccess)

                                    if (.not. tSuccess) then
                                        amp = -cc_ops(1)%get_amp(ia) * cc_ops(2)%get_amp(jk_bc) &
                                              + cc_ops(1)%get_amp(ib) * cc_ops(2)%get_amp(jk_ac) &
                                              - cc_ops(1)%get_amp(ic) * cc_ops(2)%get_amp(jk_ab) &
                                              + cc_ops(1)%get_amp(ja) * cc_ops(2)%get_amp(ik_bc) &
                                              - cc_ops(1)%get_amp(jb) * cc_ops(2)%get_amp(ik_ac) &
                                              + cc_ops(1)%get_amp(jc) * cc_ops(2)%get_amp(ik_ab) &
                                              - cc_ops(1)%get_amp(ka) * cc_ops(2)%get_amp(ij_bc) &
                                              + cc_ops(1)%get_amp(kb) * cc_ops(2)%get_amp(ij_ac) &
                                              - cc_ops(1)%get_amp(kc) * cc_ops(2)%get_amp(ij_ab) &
                                              - cc_ops(1)%get_amp(ia) * cc_ops(1)%get_amp(jb) * cc_ops(1)%get_amp(kc) &
                                              + cc_ops(1)%get_amp(ia) * cc_ops(1)%get_amp(jc) * cc_ops(1)%get_amp(kb) &
                                              + cc_ops(1)%get_amp(ib) * cc_ops(1)%get_amp(ja) * cc_ops(1)%get_amp(kc) &
                                              - cc_ops(1)%get_amp(ib) * cc_ops(1)%get_amp(jc) * cc_ops(1)%get_amp(ka) &
                                              - cc_ops(1)%get_amp(ic) * cc_ops(1)%get_amp(ja) * cc_ops(1)%get_amp(kb) &
                                              + cc_ops(1)%get_amp(ic) * cc_ops(1)%get_amp(jb) * cc_ops(1)%get_amp(ka)

                                        if (abs(amp) > EPS) then
                                            ! then add this to the list!
                                            temp_amps(n) = amp
                                            temp_ops(n, 1, :) = [elec_i, elec_j, elec_k]
                                            temp_ops(n, 2, :) = [orb_a, orb_b, orb_c]
                                            n = n + 1

                                            cc_amp_norm(0, 3) = cc_amp_norm(0, 3) + 1
                                            cc_amp_norm(1, 3) = cc_amp_norm(1, 3) + abs(amp)
                                            cc_amp_norm(2, 3) = cc_amp_norm(2, 3) + amp**2
                                        end if
                                    end if
                                end do
                            end do
                        end do
                    end do
                end do
            end do
!             end if

            cc_amp_norm(2, 3) = sqrt(cc_amp_norm(2, 3))

            print *, "triples after T1*T2 and T1^3: ", n - 1, "on proc: ", &
                iProcIndex
            ! and then allocate the actual array:
            allocate(cc_ops(3)%amplitudes(n - 1))
            allocate(cc_ops(3)%operators(n - 1, 2, 3))

            cc_ops(3)%amplitudes = temp_amps(1:n - 1)
            cc_ops(3)%operators = temp_ops(1:n - 1, :, :)
            cc_ops(3)%n_ops = n - 1
        end if

        if (cc_order >= 4) then
            ! to compare the fciqmc amplitudes and the CC amplitudes I
            ! also need to store the quadrupels of the wavefunction
            ! NOTE: for now only do that if there are no singles!
            ! as this allows us to only consider the T_2^2 influence!
            if (n_singles /= 0) then
                call stop_all(this_routine, &
                              "T_4 only implemented for zero singles!")
            end if

            if (t_store_hash_quadrupels) then
                quad_hash_size_wf = n_doubles * (n_doubles - 1) / 2
                call init_cc_hash(quad_hash_wf, quad_hash_size_wf)
            end if

            do idet = 1, int(totwalkers)

                ic = FindBitExcitLevel(iLutRef(:, 1), CurrentDets(:, idet))
                call extract_sign(CurrentDets(:, idet), sign_tmp)

                if (abs(sign_tmp(1)) > EPS) then
                    if (ic == 4) then
                        call get_bit_excitmat(iLutRef(:, 1), CurrentDets(:, idet), ex, ic)

                        ! convert it to the unique quad-index
                        quad_ind = [ex(1, 1:2), ex(2, 1:2), ex(1, 3:4), ex(2, 3:4)]
                        ! for some strange reason i decided to have a different
                        ! ordering for this than above.. change that! TODO
                        temp_int = apply_excit_ops(iLutRef(:, 1), &
                                                   [ex(1, :)], [ex(2, :)])

                        ! this should never find..
                        call cc_hash_look_up(quad_ind, temp_int, quad_hash_wf, &
                                             hash_ind, t_found)

                        if (t_found) then
                            call stop_all(this_routine, &
                                          "repeated quad not possible from wavefunction!")
                        end if

                        call cc_hash_add(quad_hash_wf, hash_ind, temp_int, sign_tmp(1))

                        cc_amp_norm(0, 4) = cc_amp_norm(0, 4) + 1

                    end if
                end if
            end do
            print *, "cc-4 L0 norm of wavefunction: ", n_quads, "on proc: ", iProcIndex
            call calc_cc_quad_norm(quad_hash_wf, quad_hash_size_wf)
            if (iProcIndex == root) then
                print *, "---test quads from wf: "
                print *, "direct sampled quads: ", n_quads
                print *, "L0 Norm of quadrupels in FCIMC-wavefunction: ", &
                    cc_amp_norm(0, 4)
                print *, "L1 Norm of quadrupels in FCIMC-wavefunction: ", &
                    cc_amp_norm(1, 4)
                print *, "L2 Norm of quadrupels in FCIMC-wavefunction: ", &
                    cc_amp_norm(2, 4)
            end if
        end if

    end subroutine fill_cc_amplitudes

    function get_ex(this, ind) result(ex)
        ! this function gives the specific excitation, if present in the
        ! cc_amplitudes, which are encoded in a linear fashion
        ! with the convention that all the electron and orbital indices are
        ! always provided in an ordered(from lowest to highest) fashion
        class(cc_amplitude), intent(in) :: this
        integer, intent(in) :: ind
        integer :: ex(2, this%order)
        character(*), parameter :: this_routine = "get_ex"

        integer :: ij, ab, cum, i, j, a, b, nij

        ASSERT(ind > 0)
        select case (this%order)
        case (1)

            if (allocated(this%operators)) then
                ex = this%operators(ind, :, :)
                return
            end if

            ! fix this below if i want to actually encode that directly and
            ! not store the operators..
            ASSERT(ind <= nel * (nbasis - nel))

            ! it is a single excition so this is easy to decompose
            ! this decomposition depends on the encoding below!
            ex(2, 1) = mod(ind - 1, (nbasis - nel))
            ! lets hope the integer division is done correctly.. on all compilers
            ex(1, 1) = int((ind - ex(2, 1)) / (nbasis - nel)) + 1
            ! and modify by nel the orbital index again to get the real
            ! orbital index!
            ex(2, 1) = ex(2, 1) + nel + 1

        case (2)

            if (allocated(this%operators)) then
                ex = this%operators(ind, :, :)
                return
            end if

            ! fix this below if i want to actually encode that directly and
            ! not store the operators..
            call stop_all(this_routine, "still buggy for double excitations!")
            ! first we have to get ij, ab back:

            nij = nel * (nel - 1) / 2

            ab = mod(ind - 1, nij)
            ij = (ind - ab) / (nij) + 1

            ab = ab + 1

            ! then we have to get ij and ab from those indices in the same
            ! way
            ! but here it gets more tricky because we can just do the mod..
            ! maybe i have to store the indices in the end..
            i = 1
            cum = (nel - i)

            ! do a stupid sum..
            do while (cum < ij .and. i <= nel)
                i = i + 1
                cum = cum + (nel - i)
            end do

            j = ij + i - (cum - nel + i)

            ! and the same for ab
            a = 1
            cum = (nbasis - nel - a)
            do while (cum < ab .and. a <= (nbasis - nel))
                a = a + 1
                cum = cum + (nbasis - nel - a)
            end do

            b = ab + a - (cum - (nbasis - nel) + a)

            ex(1, :) = [i, j]
            ex(2, :) = [a, b] + nel

        case (3)
            ex = this%operators(ind, :, :)

        case default
            call stop_all(this_routine, "higher than cc_order 2 not yet implemented!")

        end select

    end function get_ex

    function get_ind(this, elec_ind, orb_ind) result(ind)
        ! depending on the cc_order this encodes all the electron and
        ! orbital indices of an excitation in a unique linear fashion
        ! we can assume the electron and orbital indices to be ordered from
        ! highest to lowest
        class(cc_amplitude), intent(in) :: this
        integer, intent(in) :: elec_ind(this%order), orb_ind(this%order)
        integer :: ind

        character(*), parameter :: this_routine = "get_ind"

        integer :: ij, ab, nij, orb(2)

        ind = -1

        ! to i want to access this with invalid excitations?
        ASSERT(all(elec_ind > 0))
        ASSERT(all(orb_ind <= nbasis))

        ! and assert ordered inputs..
        ! or do we want to order it here, if it is not ordered?
        ! tbd!
        ASSERT(minloc(elec_ind, 1) == 1)
        ASSERT(maxloc(elec_ind, 1) == this%order)

        ASSERT(minloc(orb_ind, 1) == 1)
        ASSERT(maxloc(orb_ind, 1) == this%order)

        select case (this%order)
        case (1)
            ! single excitation
            ! the elec_ind goes from 1 to nel and the orb_ind goes from
            ! nel + 1 to nbasis (has nbasis - nel) possible values
            ! can we assume a closed shell ordered reference?
            ! with the nel lowest orbitals occupied?
            ! otherwise this is a bit more tricky..
            ind = ind_matrix_singles(elec_ind(1), orb_ind(1))
            return

            ind = (elec_ind(1) - 1) * (nbasis - nel) + (orb_ind(1) - nel)

        case (2)

            ij = linear_elec_ind(elec_ind(1), elec_ind(2))
            ab = linear_orb_ind(orb_ind(1), orb_ind(2))

            if (ij == 0 .or. ab == 0) then
                ind = 0
            else
                ind = ind_matrix_doubles(ij, ab)
            end if

            return

            call stop_all(this_routine, "still buggy for double excitations!")
            ! double excitation
            ! first encode the ij electron index in a linear fashion
            ASSERT(elec_ind(1) < elec_ind(2))
            ASSERT(orb_ind(1) < orb_ind(2))

            orb = orb_ind - nel

            ij = (elec_ind(1) - 1) * nel - elec_ind(1) * (elec_ind(1) - 1) / 2 + (elec_ind(2) - elec_ind(1))
!             ij = (elec_ind(2) - 1)*elec_ind(2)/2 + elec_ind(1) - 1
            ! i shift the orb_indices by nel..
!             ab = (orb_ind(1) - nel - 1)*(nbasis - nel) + (orb_ind(2) - orb_ind(1))
            ab = (orb(1) - 1) * (nbasis - nel) - orb(1) * (orb(1) - 1) / 2 + (orb(2) - orb(1))
!             ab = (orb_ind(2) - nel - 1)*(orb_ind(2)-nel)/2 + orb_ind(1) - nel -1

            nij = nel * (nel - 1) / 2
            ind = (ij - 1) * nij + ab

        case default
            call stop_all(this_routine, "higher than cc_order 2 not yet implemented!")
        end select

    end function get_ind


    pure integer function ind1(i, a)

        implicit none
        integer, intent(in) :: i, a

        ind1 = (i - 1) * (nbasis - nel) + a

        return
    end function ind1

    pure integer function ind2(i, j, a, b)

        implicit none
        integer, intent(in) :: i, j, a, b
        integer :: ij, ab, nvirt, nab

        ij = (j - 1) * j / 2 + i
        ab = (b - 1) * b / 2 + a
        nvirt = nbasis - nel
        nab = nvirt * (nvirt + 1) / 2

        ind2 = (ij - 1) * nab + ab

        return
    end function ind2

    function calc_number_of_excitations(n_alpha, n_beta, max_excit, n_orbs) &
        result(n_excits)
        ! i might want a routine which calculates the number of all possible
        ! excitations for each excitation level
        integer, intent(in) :: n_alpha, n_beta, max_excit, n_orbs
        integer :: n_excits(max_excit)

        integer :: n_parallel(max_excit, 2), i, j, k

        ! the number of all possible single excitations, per spin is:

        n_parallel = 0
        n_excits = 0

        do i = 1, max_excit
            n_parallel(i, 1) = calc_n_parallel_excitations(n_alpha, n_orbs, i)
            n_parallel(i, 2) = calc_n_parallel_excitations(n_beta, n_orbs, i)
        end do

        ! i can set this up in a general way..

        do i = 1, max_excit

            ! the 2 parallel spin species always are added
            n_excits(i) = sum(n_parallel(i, :))

            ! and i have to zig-zag loop over the other possible combinations
            ! to lead to this level of excitation..
            j = i - 1
            k = 1

            do while (j >= k)

                n_excits(i) = n_excits(i) + n_parallel(j, 1) * n_parallel(k, 2)

                ! and if j /= k do it the other way around too
                if (j /= k) then
                    n_excits(i) = n_excits(i) + n_parallel(j, 2) * n_parallel(k, 1)

                end if
                ! and in/decrease the counters
                j = j - 1
                k = k + 1

                ! in this way i calculate eg.:
                ! n_ij^ab = n_ij_ab(u + d) + n_i^a(u) * n_i^a(d)
                ! and
                ! n_ijk^abc = n_ijk^abc(u + d) + n_ij^ab(u) * n_i^a(d) and spin-flipped..

            end do
        end do

    end function calc_number_of_excitations

    function calc_n_parallel_excitations(n_elecs, n_orbs, ic) result(n_parallel)
        ! this function determines the number of parallel spin excitaiton
        ! with each electrons having the same spin for a given excitation
        ! level and number of electrons and available orbitals for this spin
        integer, intent(in) :: n_elecs, n_orbs, ic
        integer :: n_parallel

        n_parallel = int(choose_i64(n_elecs, ic) * choose_i64(n_orbs - n_elecs, ic))

    end function calc_n_parallel_excitations
end module cc_amplitudes