CountExcitations3 Subroutine

    SUBROUTINE CountExcitations3(nI, exflag, nSingleExcits, nDoubleExcits)
! This routine simply counts the excitations in terms of single and doubles from the nI determinant.
! The exflag sent through indicates which should be counted - exflag=1 means only singles, exflag=2 means
! only doubles, and anything else both are counted.
        USE SymData, only: nSymLabels
        USE SystemData, only: ElecPairs, tFixLz, iMaxLz
        USE GenRandSymExcitNUMod, only: PickElecPair, construct_class_counts, ClassCountInd, ScratchSize
        INTEGER :: nSingleExcits, nDoubleExcits, Symi, i, Spini, nI(NEl)
        INTEGER :: iSpn, Elec1Ind, Elec2Ind, SymProduct, exflag
        INTEGER :: Syma, Symb, Spina, Spinb, StartSpin, EndSpin
        INTEGER :: ClassCount2(ScratchSize), SumMl
        INTEGER :: ClassCountUnocc2(ScratchSize)
        INTEGER :: StartMl, EndMl, Mla, Mlb

        CALL construct_class_counts(nI, ClassCount2, ClassCountUnocc2)
! This sets up arrays containing the number of occupied and unoccupied in each symmetry.
! ClassCounts2(1,:)=No alpha occupied, ClassCounts2(2,:)=No Beta occupied.
! ClassCountsUnocc2(1,:)=No alpha unocc, ClassCounts2Unocc2(2,:)=No Beta unocc.
! The second index of these arrays referrs to the symmetry (0 -> 7).

! Only counting.  Run through each occupied orbital, and count the number of spin and symmetry allowed orbitals it
! may be excited to.
        nSingleExcits = 0
        nDoubleExcits = 0

        IF (exflag /= 2) THEN
! Count the singles.
! Take each electron and find out the number of symmetry allowed orbitals it may be excited to.
            do i = 1, NEl
                Symi = SpinOrbSymLabel(nI(i))
                IF ((G1(nI(i))%Ms) == -1) Spini = 2        ! G1(i)%Ms is -1 for beta, and 1 for alpha.
                IF ((G1(nI(i))%Ms) == 1) Spini = 1         ! Translate this into 1 for alpha and 2 for beta
                ! for the ClassCount arrays.
                IF (tFixLz) THEN
                    Mla = G1(nI(i))%Ml
                ELSE
                    Mla = 0
                end if

! This electron in orbital of SymI and SpinI can only be excited to orbitals with the same spin and symmetry.
! Then add in the number of unoccupied orbitals with the same spin and symmetry to which each electron may be excited.

                nSingleExcits = nSingleExcits + ClassCountUnocc2(ClassCountInd(Spini, Symi, Mla))

            end do
        end if

! This is the end of the singles.
!            write(stdout,*) 'Number of singles',nSingleExcits

! For the doubles, first pick an electron pair i,j.
! Based on these orbitals, run through each spin and each symmetry - take this to be orbital a.
! Multiply the number with these symmetries by the number of possible b orbitals which correspond.
! Do this for all a and then all i,j pairs.

        IF (exflag /= 1) THEN
            do i = 1, ElecPairs

! iSpn=2 for alpha beta pair, ispn=3 for alpha alpha pair and ispn=1 for beta beta pair.
                CALL PickElecPair(nI, Elec1Ind, Elec2Ind, SymProduct, iSpn, SumMl, i)

                StartSpin = 1
                EndSpin = 2
                IF (iSpn == 3) EndSpin = 1
                IF (iSpn == 1) StartSpin = 2
                do Spina = StartSpin, EndSpin            ! Run through both spins, orbital a may be alpha or beta.
                    IF (iSpn == 2) THEN
! Spin of orbital b should be opposite to orbital a.
                        IF (Spina == 1) Spinb = 2
                        IF (Spina == 2) Spinb = 1
                    ELSE
! Spin of orbital b should be the same as orbital a.
                        IF (Spina == 1) Spinb = 1
                        IF (Spina == 2) Spinb = 2
                    end if

                    do Syma = 0, nSymLabels - 1

! Need to work out the symmetry of b, given the symmetry of a (Sym).
                        Symb = IEOR(Syma, SymProduct)

                        IF (tFixLz) THEN
                            StartMl = -iMaxLz
                            EndMl = iMaxLz
                        ELSE
                            StartMl = 0
                            EndMl = 0
                        end if

                        do Mla = StartMl, EndMl

                            Mlb = SumMl - Mla   !Will be 0 if no Lz, otherwise we need Mla + Mlb = Mli + Mlj = SumMl

                            IF (ABS(Mlb) <= iMaxLz) THEN
                                IF ((Spina == Spinb) .and. (Syma == Symb) .and. (Mla == Mlb)) THEN
                                    ! If the spin and spatial symmetries of a and b are the same
                                    ! there will exist a case where Orba = Orbb, want to remove this.
                                    nDoubleExcits = nDoubleExcits + (ClassCountUnocc2(ClassCountInd(Spina, Syma, Mla)) &
                                                                     * (ClassCountUnocc2(ClassCountInd(Spinb, Symb, Mlb)) - 1))
                                ELSE
                                    nDoubleExcits = nDoubleExcits + (ClassCountUnocc2(ClassCountInd(Spina, Syma, Mla)) &
                                                                     * ClassCountUnocc2(ClassCountInd(Spinb, Symb, Mlb)))
                                end if
                            end if
                        end do

                    end do

                end do
            end do
            nDoubleExcits = nDoubleExcits / 2

        end if

    ENDSUBROUTINE CountExcitations3