CalcAllab Subroutine

    SUBROUTINE CalcAllab(nI, ILUT, Elec1Ind, Elec2Ind, SymProduct, iSpn, OrbA, OrbB, nParts, iCreate, Tau)
        INTEGER :: nI(NEl), Elec1Ind, Elec2Ind, SymProduct, iSpn, OrbA, OrbB, iCreate
        INTEGER(KIND=n_int) :: iLut(0:NIfTot)
        INTEGER :: SpatOrbi, SpatOrbj, Spini, Spinj, i, aspn, bspn, SymA, SymB, SpatOrba, EndSymState, VecInd
        real(dp) :: Tau, SpawnProb(MaxABPairs), NormProb, rat, r
        INTEGER :: SpawnOrbs(2, MaxABPairs), j, nParts, SpinIndex, Ind
        HElement_t(dp) :: HEl

!We want the spatial orbital number for the ij pair (Elec1Ind is the index in nI).
!Later, we'll have to use GTID for UHF.
        SpatOrbi = ((nI(Elec1Ind) - 1) / 2) + 1
        SpatOrbj = ((nI(Elec2Ind) - 1) / 2) + 1
        Spini = G1(nI(Elec1Ind))%Ms
        Spinj = G1(nI(Elec2Ind))%Ms
        VecInd = 1
        NormProb = 0.0_dp

        do i = 1, nBasis
!Run through all a orbitals
            IF (mod(i, 2) == 0) THEN
!We have an alpha spin...
                aspn = 1
                IF (iSpn == 1) THEN
!ij is an beta/beta pair, so we only want to pick beta a orbitals.
                    CYCLE
                end if
            ELSE
!We have a beta spin...
                aspn = -1
                IF (iSpn == 3) THEN
!ij is a alpha/alpha pair, so we only want to pick alpha a orbitals.
                    CYCLE
                end if
            end if

!We also want to check that the a orbital we have picked is not already in the determinant we are exciting from.
            IF (BTEST(ILUT((i - 1) / bits_n_int), MOD((i - 1), bits_n_int))) THEN
!Orbital is in nI...do not make a pair from this.
                CYCLE
            end if

!Now we have to run over all b's which can go with a. We only want unique pairs, so we constrain b to be at least a+1.
!We only want to run over symmetry and spin allowed b's though.
!Find the required symmetry of b.
            SymA = INT(G1(i)%Sym%S, 4)
            SymB = IEOR(SymA, SymProduct)
            SpatOrba = ((i - 1) / 2) + 1

!We also want to take into account spin.
            IF (ispn == 1) THEN
                bspn = -1  !Want beta spin b orbitals
                SpinIndex = 2
            else if (ispn == 3) THEN
                bspn = 1  !Want alpha spin b orbitals
                SpinIndex = 1
            ELSE
!ij pair is an alpha/beta spin pair, therefore b wants to be of opposite spin to a.
                IF (aspn == -1) THEN
!a orbital is a beta orbital, therefore we want b to be an alpha orbital.
                    bspn = 1
                    SpinIndex = 1
                ELSE
                    bspn = -1
                    SpinIndex = 2
                end if
            end if

!To run just through the states of the required symmetry we want to use SymLabelCounts.
!            StartSymState=SymLabelCounts(1,SymB+1)
            Ind = ClassCountInd(SpinIndex, SymB, 0)
            EndSymState = SymLabelCounts2(1, Ind) + SymLabelCounts2(2, Ind) - 1

!Run over all possible b orbitals
            do j = SymLabelCounts2(1, Ind), EndSymState

!                IF(bspn.eq.-1) THEN
!                    OrbB=(2*SymLabelList(j))-1     !This is the spin orbital chosen for b
!                ELSE
!                    OrbB=(2*SymLabelList(j))
!                end if

                OrbB = SymLabelList2(j)     !This is the spin orbital chosen for b

                IF (OrbB <= i) THEN
!Since we only want unique ab pairs, ensure that b > a.
                    CYCLE
                end if

                IF (BTEST(ILUT((OrbB - 1) / bits_n_int), MOD((OrbB - 1), bits_n_int))) THEN
!Orbital is in nI...do not make a pair from this.
                    CYCLE
                end if

!We have now found an allowed ab pair to go with the ij pair chosen previously - record its stats.
                IF (Spini == aspn .and. Spinj == bspn) THEN
                    Hel = get_umat_el(SpatOrbi, SpatOrbj, SpatOrba, j)
                ELSE
                    Hel = (0.0_dp)
                end if
                IF (Spini == bspn .and. Spinj == aspn) THEN
                    Hel = Hel - get_umat_el(SpatOrbi, SpatOrbj, j, SpatOrba)
                end if

                SpawnProb(VecInd) = abs(REAL(Hel, dp))
                SpawnOrbs(1, VecInd) = i
                SpawnOrbs(2, VecInd) = OrbB
                NormProb = NormProb + SpawnProb(VecInd)
                VecInd = VecInd + 1
                IF (VecInd >= MaxABPairs) THEN
                    CALL Stop_All("CalcAllab", "Finding too many ab pairs...")
                end if

            end do
        end do

        VecInd = VecInd - 1     !This now indicates the total number of ab pairs we have found for the chosen ij.

        IF (VecInd == 0) THEN
            CALL Stop_All("CalcAllab", "No ab pairs found for the chosen ij")
        end if

!We now have to find out how many children to spawn, based on the value of normprob.
        rat = Tau * NormProb * REAL(ElecPairs * nParts, dp) / PDoubNew
        iCreate = INT(rat)
        rat = rat - REAL(iCreate, dp)
        r = genrand_real2_dSFMT()
        IF (rat > r) THEN
!Child is created
            iCreate = iCreate + 1
        end if

        IF (iCreate > 0) THEN
!We want to spawn particles. This only question now is where. Run through the
!ab pairs again and choose based on the SpawnProb element.
            r = genrand_real2_dSFMT()
            r = r * NormProb

            i = 0
            do while (r > 0.0_dp)
                i = i + 1
                r = r - SpawnProb(i)
            end do
            IF (i > VecInd) THEN
                CALL Stop_All("CalcAllab", "Chosen ab pair does not correspond to allowed pair")
            end if

            OrbA = SpawnOrbs(1, i)
            OrbB = SpawnOrbs(2, i)

        end if

    END SUBROUTINE CalcAllab