assign_reference_dets Subroutine

    subroutine assign_reference_dets()

        ! Depending on the configuration we may have one, or multiple,
        ! reference determinants.

        integer :: det(nel), orbs(nel), orb, orb2, norb
        integer :: run, cc_idx, label_idx, i, j, found_orbs(inum_runs)
        real(dp) :: energies(nel), hdiag

        ! for now guga only works with non-complex code
        integer(n_int), allocatable :: excitations(:, :)
        integer :: n_excits, ierr
        real(dp), allocatable :: diag_energies(:)
        logical, allocatable :: found_mask(:)
        character(*), parameter :: this_routine = "assign_reference_dets"

        ! If the user has specified all of the (multiple) reference states,
        ! then just copy these across to the ilutRef array:
        if (tMultipleInitialStates) then

            tReplicaReferencesDiffer = .true.

            do run = 1, inum_runs
                ProjEDet(:, run) = initial_states(:, run)
                call EncodeBitDet(ProjEDet(:, run), ilutRef(:, run))
            end do

        else if (tOrthogonaliseReplicas) then

            tReplicaReferencesDiffer = .true.

            ! The first replica is just a normal FCIQMC simulation.
            ilutRef(:, 1) = ilutHF
            ProjEDet(:, 1) = HFDet

            found_orbs = 0

            ! in the GUGA approach we have to change that simple
            ! single excitation of course or otherwise we get non-
            ! allowed CSF or the wrong STOT symmetry sector.
            if (tGUGA) then
                ! run the exact single excitations on the HF det
                ! and find inum_runs lowest energetically ones..

                ! create all excitations from the HF
                ! do i remain in the same symmetry sector with this
                ! info??
                ! i think my exact hamil application routine does not
                ! deal with symmetry at all...
                ! i guess i have to change that for the real-space
                ! hubbard model implementation!
                if (tHUB .or. tUEG .or. .not. (tNoBrillouin)) then
                    call actHamiltonian(ilutHF, CSF_Info_t(ilutHF), excitations, n_excits)
                else
                    call actHamiltonian(ilutHF, CSF_Info_t(ilutHF), excitations, n_excits, .true.)
                end if

                ! if no excitations possible... there is something wrong
                if (.not. (n_excits > 0) .or. n_excits < inum_runs - 1) then
                    print *, "n_excits: ", n_excits
                    print *, "requested excited states:", inum_runs - 1
                    call write_guga_list(stdout, excitations(:, 1:n_excits))
                    call stop_all(this_routine, "not enough excitations from HF!")
                end if

                ! and the choose the inum_runs - 1 lowest energetically
                ! excitations

                ! also have to calculate the correct diagonal element to
                ! compare energies
                allocate(diag_energies(n_excits), stat=ierr)
                allocate(found_mask(n_excits), stat=ierr)
                found_mask = .true.

                do i = 1, n_excits
                    diag_energies(i) = real(calcDiagMatEleGUGA_ilut(excitations(:, i)), dp)
                end do

                ! can i sort the excitation list, according to energies?
                do run = 2, inum_runs
                    ! find the minimum energy
                    i = minloc(diag_energies, 1, found_mask)
                    ! dont find that specific one again
                    found_mask(i) = .false.
                    ! assign the reference determinant(transform between
                    ! guga and actual iluts!)
                    call convert_ilut_toNECI(excitations(:, i), ilutRef(:, run))
                    ! and calc. the nI representation
                    call decode_bit_det(ProjEDet(:, run), ilutRef(:, run))
                    ! that should be it.. should the diagonal element also
                    ! be already stored within the ilut? is that done here?
                end do

            else
                do run = 2, inum_runs

                    ! Now we want to find the lowest energy single excitation with
                    ! the same symmetry as the reference site.

                    do i = 1, nel
                        ! Find the excitations, and their energy
                        orb = HFDet(i)
                        cc_idx = ClassCountInd(orb)
                        label_idx = SymLabelCounts2(1, cc_idx)
                        norb = OrbClassCount(cc_idx)

                        ! nb. sltcnd_0 does not depend on the ordering of the det,
                        !     so we don't need to do any sorting here.
                        energies(i) = 9999999.9_dp
                        do j = 1, norb
                            orb2 = SymLabelList2(label_idx + j - 1)
                            if (.not. (any(orb2 == HFDet) .or. any(orb2 == found_orbs))) then
                                det = HFDet
                                det(i) = orb2
                                hdiag = real(sltcnd_excit(det, Excite_0_t()), dp)
                                if (hdiag < energies(i)) then
                                    energies(i) = hdiag
                                    orbs(i) = orb2
                                end if
                            end if
                        end do
                    end do

                    ! Which of the electrons that is excited gives the lowest energy?
                    i = minloc(energies, 1)
                    found_orbs(run) = orbs(i)

                    ! Construct that determinant, and set it as the reference.
                    ProjEDet(:, run) = HFDet
                    ProjEDet(i, run) = orbs(i)
                    call sort(ProjEDet(:, run))
                    call EncodeBitDet(ProjEDet(:, run), ilutRef(:, run))

                end do
            end if

        else if (tPreCond) then

            do run = 1, inum_runs
                ilutRef(:, run) = ilutHF
                ProjEDet(:, run) = HFDet
            end do

            ! And make sure that the rest of the code knows this
            tReplicaReferencesDiffer = .true.

        else

            ! This is the normal case. All simultions are essentially doing
            ! the same thing...

            do run = 1, inum_runs
                ilutRef(:, run) = ilutHF
                ProjEDet(:, run) = HFDet
            end do

            ! And make sure that the rest of the code knows this
            tReplicaReferencesDiffer = .false.

        end if

        write(stdout, *) 'Generated reference determinants:'
        do run = 1, inum_runs
            call write_det(stdout, ProjEDet(:, run), .false.)
            hdiag = real(get_helement(ProjEDet(:, run), ProjEDet(:, run), 0), dp)
            write(stdout, '(" E = ", f16.9)') hdiag
        end do

    end subroutine assign_reference_dets