DetInit Module Subroutine

    module subroutine DetInit()
        real(dp) :: DNDET
        integer :: i, j
        integer(int64) :: nDet
        integer :: alpha, beta, symalpha, symbeta, endsymstate
        LOGICAL :: tSuccess, tFoundOrbs(nBasis)

        write (stdout, *) "SYMMETRY MULTIPLICATION TABLE"
        CALL WRITESYMTABLE(stdout)

        CALL GENSymStatePairs(NBASIS / 2, .false.)

!iActiveBasis is a copy of nPaths
        IF (iActiveBasis == -2) then
!  PATHS ACTIVE SETS
            Call GenActiveBasis(ARR, nBasis, nEl, nActiveBasis, nActiveSpace(1), nActiveSpace(2))
        else if (iActiveBasis == -3) then
!  PATHS ACTIVE ORBITALS
            nActiveBasis(1) = nEl + 1 - nActiveSpace(1)
            nActiveBasis(2) = nEl + nActiveSpace(2)
            write (stdout, *) "Active space:", nActiveBasis(1), " TO ", nActiveBasis(2), " (ordered labels)."
            write (stdout, *) "Active space electrons:", nEl - nActiveBasis(1) + 1
        else
            nActiveBasis(1) = 1
            nActiveBasis(2) = nBasis
        end if
!C.. Work out a preliminary Fermi det
!      IF(FDET(1).EQ.0) THEN

!C.. Check if we're blocking the hamiltonian
!C      IF(THFBASIS.AND.TBLOCK) THEN
!C         write(stdout,*) "THFBASIS set and NBLK=0.  ",
!C     &         "Cannot block diagonalize in HF Basis."
!C         STOP
!C      end if
!C      CALL SYMGENEXCITS(FDET,NEL,NBASIS)
!C      CALL LeaveMemoryManager
!C      STOP

!C.. in order to calculate the H matrix, we need to work out all the determinants
!C.. beware with NPATHS - it calcs the list of dets even if we don't calc H
!C.. Could be big.
!C..Now we see how many determinants we need
!C      IF(nBasis.GT.170) THEN
!C..This fix is to stop floating overflow as taking the factorial of (nBasis.GT.170) crashes
!C  using the old FACTRL routine.
        NDET = 1
        DNDET = 1.0_dp
        DO I = 0, NEL - 1
            NDET = (NDET * (nBasis - I)) / (I + 1)
            DNDET = (DNDET * real(nBasis - I, dp)) / real(I + 1, dp)
        end do

        IF (abs(real(NDET) - dndet) > 1.0e-6) THEN
!         write(stdout,*) ' NUMBER OF DETERMINANTS : ' , DNDET
            NDET = -1
        ELSE
!         write(stdout,*) ' NUMBER OF DETERMINANTS : ' , NDET
        end if

!C      CALL TC(I_HMAX,I_P,NWHTAY)

!Check that the symmetry routines have set the symmetry up correctly...
        tSuccess = .true.
        tFoundOrbs(:) = .false.

        IF ((.not. tHub) .and. (.not. tUEG) .and. TwoCycleSymGens) THEN
            do i = 1, nSymLabels
                EndSymState = SymLabelCounts(1, i) + SymLabelCounts(2, i) - 1
                do j = SymLabelCounts(1, i), EndSymState

                    Beta = (2 * SymLabelList(j)) - 1
                    Alpha = (2 * SymLabelList(j))
                    SymAlpha = INT((G1(Alpha)%Sym%S), 4)
                    SymBeta = INT((G1(Beta)%Sym%S), 4)

                    IF (.not. tFoundOrbs(Beta)) THEN
                        tFoundOrbs(Beta) = .true.
                    ELSE
                        CALL Stop_All("SetupParameters", "Orbital specified twice")
                    end if
                    IF (.not. tFoundOrbs(Alpha)) THEN
                        tFoundOrbs(Alpha) = .true.
                    ELSE
                        CALL Stop_All("SetupParameters", "Orbital specified twice")
                    end if

                    IF (G1(Beta)%Ms /= -1) THEN
                        tSuccess = .false.
                    else if (G1(Alpha)%Ms /= 1) THEN
                        tSuccess = .false.
                    else if ((SymAlpha /= (i - 1)) .or. (SymBeta /= (i - 1))) THEN
                        tSuccess = .false.
                    end if
                end do
            end do
            do i = 1, nBasis
                IF (.not. tFoundOrbs(i)) THEN
                    write (stdout, *) "Orbital: ", i, " not found."
                    CALL Stop_All("SetupParameters", "Orbital not found")
                end if
            end do
        end if
        ! SpinOrbSymSetup currently sets up the symmetry arrays for use with
        ! symrandexcit2 excitation routines. These are not currently
        ! compatible with non-abelian symmetry groups, which CPMD jobs
        ! invariably used. To avoid this complication, this symmetry
        ! setup will not be used with CPMD, and thus these excitation
        ! generators won't work.
        IF (.not. tCPMD) THEN
            IF (.not. tSuccess) THEN
                write (stdout, *) "************************************************"
                write (stdout, *) "**                 WARNING!!!                 **"
                write (stdout, *) "************************************************"
                write (stdout, *) "Symmetry information not set up correctly in NECI initialisation"
                write (stdout, *) "Will attempt to set up the symmetry again, but now in terms of spin orbitals"
                write (stdout, *) "Old excitation generators will not work"
                write (stdout, *) "I strongly suggest you check that the reference energy is correct."
                !CALL SpinOrbSymSetup() !.true.)
            ELSE
                write (stdout, *) "Symmetry and spin of orbitals correctly set up for excitation generators."
                write (stdout, *) "Simply transferring this into a spin orbital representation."
                !CALL SpinOrbSymSetup() !.false.)
            end if

            if (tPickVirtUniform) then
                call virt_uniform_sym_setup()
            else
                ! Includes normal & HPHF
                call SpinOrbSymSetup()
            end if

        end if
! From now on, the orbitals are also contained in symlabellist2 and symlabelcounts2.
! These are stored using spin orbitals.
    end subroutine DetInit