DetPreFreezeInit Subroutine

    Subroutine DetPreFreezeInit()

        integer ierr, ms, iEl, flagAlpha
        integer i, j, Lz, OrbOrder(8, 2), FDetTemp(NEl), lmsMax
        type(BasisFn) s
        logical :: tGenFDet
        HElement_t(dp) :: OrbE
        character(25), parameter :: this_routine = 'DetPreFreezeInit'
        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

            ! A quick check that we have defined a reasonable det.
            ms = sum(get_spin_pn(fdet(1:nel)))
            tRef_Not_HF = .true.
        else
            tGenFDet = .true.
            lmsMax = sum(nOccOrbs) - sum(nClosedOrbs)
            if ((sum(nOccOrbs) + sum(nClosedOrbs)) == nel &
                .and. (.not. TSPN .or. abs(LMS) == lmsMax)) then
                tGenFDet = .false.
                if (LMS < 0) then
                    flagAlpha = 0
                else
                    flagAlpha = 1
                end if
                iEl = 1
                do i = 1, nIrreps
                    ! nClosedOrbs is the number of alpha/beta orbs (the ones with minority spin)
                    do j = 1, nClosedOrbs(i)
                        FDet(iEl) = irrepOrbOffset(i) + 2 * j - flagAlpha
                        iEl = iEl + 1
                    end do
                    ! nOccOrbs is the number of majority spin orbs (per irrep)
                    do j = 1, nOccOrbs(i)
                        FDet(iEl) = irrepOrbOffset(i) + 2 * j - (1 - flagAlpha)
                        iEl = iEl + 1
                    end do
                end do
                call sort(FDet)
            end if
            if (tGenFDet) then
                CALL GENFDET(FDET)
                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 .or. tZeroDelimitedFCIDUMP) then
            !Orbitals are ordered by occupation number from MOLPRO, and not reordered in NECI
            !Therefore, we know HF determinant is first four occupied orbitals.
            if (tMolpro) 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
            if (tMolpro .or. tHFOrder .or. .not. tHFNoOrder) then
                write (stdout, "(A)") "Reordering basis by orbital energies..."
                OrbOrder(:, :) = 0
                call ORDERBASIS(NBASIS, Arr, Brr, OrbOrder, nBasisMax, G1)
            else
                write (stdout, "(A)") "Not reordering basis."
                Arr(:nBasis, 2) = Arr(:nBasis, 1)
            end if
            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)))

    End Subroutine DetPreFreezeInit