DetPreFreezeInit_old Subroutine

    Subroutine DetPreFreezeInit_old()
        integer :: ierr, ms, iEl, flagAlpha, msTmp
        integer :: i, j, Lz, OrbOrder(8, 2), FDetTemp(NEl)
        type(BasisFn) s
        logical :: tGenFDet
        HElement_t(dp) :: OrbE
        character(25), parameter :: this_routine = 'DetPreFreezeInit_old'
        allocate (FDet(nEl), stat=ierr)
        LogAlloc(ierr, 'FDet', nEl, 4, tagFDet)
        IF (tDefineDet) THEN
            write (stdout, *) 'Defining FDet according to input'
            do i = 1, NEl
                FDet(i) = DefDet(i)
            end do
            call assignOccOrbs()

            ! A quick check that we have defined a reasonable det.
            ms = sum(get_spin_pn(fdet(1:nel)))
            if (abs(ms) /= abs(lms)) then
                write (stdout, *) 'LMS', lms
                write (stdout, *) 'Calculated Ms', ms
                call stop_all(this_routine, "Defined determinant has the &
                     &wrong Ms value. Change DEFINEDET or &
                     &SPIN-RESTRICT")
            end if

            if (tGUGA) then
                ms = abs(ms)
                if (ms < 0 .or. ms /= STOT) then
                    write (stdout, *) "S: ", STOT
                    write (stdout, *) "calculated S of inputted CSF: ", ms
                    call stop_all(this_routine, " Defined CSF has the &
                        &wrong total spin quantum number. Change DEFINEDET or &
                        &S quantum numnber!")
                end if
            end if

            tRef_Not_HF = .true.
        else
            if ((sum(nOccOrbs) + sum(nClosedOrbs)) == nel) then
                tGenFDet = .false.
                iEl = 1
                msTmp = -1 * lms
                do i = 1, nIrreps
                    ! doubly occupy the closed orbs
                    do j = 1, nClosedOrbs(i)
                        FDet(iEl) = irrepOrbOffset(i) + 2 * j - 1
                        iEl = iEl + 1
                        FDet(iEl) = irrepOrbOffset(i) + 2 * j
                        iEl = iEl + 1
                    end do
                    ! now distribute electrons to the open orbs
                    ! such that the total sz matches the requested
                    ! only consider occ orbs which are not closed
                    do j = nClosedOrbs(i) + 1, nOccOrbs(i)
                        if (msTmp < 0) then
                            flagAlpha = 0
                        else
                            flagAlpha = 1
                        end if
                        FDet(iEl) = irrepOrbOffset(i) + 2 * j - flagAlpha
                        iEl = iEl + 1
                        msTmp = msTmp + (1 - 2 * flagAlpha)
                    end do
                end do
                call sort(FDet)
            else
                CALL GENFDET(FDET)
                call assignOccOrbs()
                IF (tUEGSpecifyMomentum) THEN
                    write (stdout, *) 'Defining FDet according to a momentum input'
                    CALL ModifyMomentum(FDET)
                end if
                tRef_Not_HF = .false.
            end if
        end if
        write (stdout, "(A)", advance='no') " Fermi det (D0):"
        call write_det(stdout, FDET, .true.)
        Call GetSym(FDet, nEl, G1, nBasisMax, s)
        write (stdout, "(A)", advance='no') " Symmetry: "
        Call writesym(stdout, s%Sym, .true.)
        IF (tFixLz) THEN
            Call GetLz(FDet, nEl, Lz)
            write (stdout, "(A,I5)") "Lz of Fermi det:", Lz
        end if
        CALL NECI_ICOPY(NEL, FDET, 1, NUHFDET, 1)
        if (tMolpro) then
            !Orbitals are ordered by occupation number from MOLPRO, and not reordered in NECI
            !Therefore, we know HF determinant is first four occupied orbitals.
            write (stdout, "(A)") "NECI called from MOLPRO, so assuming orbitals ordered by occupation number."
            if (.not. tDefineDet) then
                FDetTemp(:) = FDet(:)
            else
                !We have defined our own reference determinant, but still use the first orbitals for the calculation
                !of 'orbital energies'
                CALL GENFDET(FDETTEMP)
            end if
            write (stdout, "(A)") "Calculating orbital energies..."
            do i = 1, nBasis
                OrbE = CalcFockOrbEnergy(i, FDetTemp)
                Arr(i, 1) = real(OrbE, dp)
                Brr(i) = i
            end do
            write (stdout, "(A)") "Reordering basis by orbital energies..."
            OrbOrder(:, :) = 0
            call ORDERBASIS(NBASIS, Arr, Brr, OrbOrder, nBasisMax, G1)
            !However, we reorder them here
            call writebasis(stdout, G1, nBasis, Arr, Brr)
        end if
        E0HFDET = ECORE
        DO I = 1, NEL
            E0HFDET = E0HFDET + ARR(NUHFDET(i), 2)
        end do
        write (stdout, *) "Fock operator energy:", E0HFDET

        ! Store the value of Ms for use in other areas
        calculated_ms = sum(get_spin_pn(fdet(1:nel)))

        if (tGUGA) then
            calculated_ms = abs(calculated_ms)
        end if

    contains

        subroutine assignOccOrbs
            integer :: k
            ! assign the occ/closed orbs
            nOccOrbs = 0
            nClosedOrbs = 0
            do k = 1, nel
                if (k > 1) then
                    if (is_alpha(FDet(k)) .and. FDet(k - 1) == FDet(k) - 1) then
                        ! we do not need to resolve the irrep anymore (not relevant
                        ! for further usage)
                        nClosedOrbs(1) = nClosedOrbs(1) + 1
                        cycle
                    end if
                end if
                nOccOrbs(1) = nOccOrbs(1) + 1
            end do
        end subroutine assignOccOrbs
    End Subroutine DetPreFreezeInit_old