init_guga Subroutine

    subroutine init_guga()
        integer :: i
        character(*), parameter :: this_routine = "init_guga"
        ! main initialization routine

        ! this routine is called in SysInit() of System_neci.F90
        write(stdout, *) ' ************ Using the GUGA-CSF implementation **********'
        write(stdout, *) ' Restricting the total spin of the system, tGUGA : ', tGUGA
        write(stdout, '(A,I5)') '  Restricting total spin S in units of h/2 to ', STOT
        write(stdout, *) ' So eg. S = 1 corresponds to one unpaired electron '
        write(stdout, *) ' not quite sure yet how to deal with extensively used m_s quantum number..'
        write(stdout, *) ' NOTE: for now, although SPIN-RESTRICT is set off, internally m_s(LMS) '
        write(stdout, *) ' is set to STOT, to make use of reference determinant creations already implemented'
        write(stdout, *) ' Since NECI always seems to take the beta orbitals first for open shell or '
        write(stdout, *) ' spin restricted systems, associate those to positively coupled +h/2 orbitals '
        write(stdout, *) ' to always ensure a S >= 0 value!'
        write(stdout, *) ' *********************************************************'

        if (tGen_nosym_guga) then
            call Stop_All(this_routine, "'nosym-guga' option deprecated!")
        end if

        if (t_fci_pchb_excitgen) then
            call stop_all(this_routine, &
                "please specify 'guga-pchb' as excitation generator to work with GUGA!")
        end if

        ! initialize the procedure pointer arrays, needed in the matrix
        ! element calculation
        call init_guga_data_procPtrs()

        ! define global variable of spatial orbitals
        ! do that in a more general setup routine! where nBasis is defined
        ! eg
        ! i have to all this routine again from a point after freezing
        ! where the new number of NBasis is determined already..
        nSpatOrbs = nBasis .div. 2
        nSpatOrbs_i8 = int(nBasis .div. 2, int64)
        ij_pairs = fuse_symm_idx(nSpatOrbs_i8, nSpatOrbs_i8)

        ! but also have to set up the global orbitalIndex list
        orbitalIndex = [(i, i = 1, nSpatOrbs)]


        ! Store GUGA specific information about the current CSF.
        ! In principle this is redundant and could be computed from nI or ilut,
        !   but we precompute it for performance reasons.
        call new_CSF_Info_t(nSpatOrbs, current_csf_i)
        if (allocated(csf_ref)) deallocate(csf_ref)
        allocate(csf_ref(inum_runs))
        call new_CSF_Info_t(nSpatOrbs, csf_ref)

        call init_guga_class()

        ! initialize and point the excitation generator functions to the
        ! correct ones
        call init_guga_orbital_pickers()

        ! also have to set tRealCoeffs true to use it in excitation creation
        ! dont actually need realCoeffs anymore since i changed the accessing
        ! of the ilut lists used for excitation creation
        ! but probably have to set a few other things so it works with
        ! reals
        tUseRealCoeffs = .true.

        tUseFlags = .true.

        ! for now (time/iteration comparison) reasons, decide which
        ! reference energy calculation method we use
        ! use the new "direct" calculation method
        calc_off_diag_guga_ref => calc_off_diag_guga_ref_direct

        ! make checks for the RDM calculation
        if (tRDMonfly) then
            call check_rdm_guga_setup()
        end if

        ! make a unified bit rep initializer:
        call init_guga_bitrep(nifd)

        ! set some defaults for non-working things:
        t_fast_pops_core = .false.
        ss_space_in%tApproxSpace = .false.
        trial_space_in%tApproxSpace = .false.

    end subroutine init_guga