CreateDoubExcit Subroutine

    SUBROUTINE CreateDoubExcit(nI, nJ, ClassCountUnocc2, ILUT, ExcitMat, tParity, pGen)
        integer, intent(in) :: nI(nel)
        integer, intent(out) :: nJ(nel), ExcitMat(2, maxExcit)
        integer, intent(in) :: ClassCountUnocc2(ScratchSize)
        integer(n_int), intent(in) :: iLut(0:NIfTot)
        real(dp), intent(out) :: pGen
        logical, intent(out) :: tParity

        integer :: nExcitOtherWay, orbB, nExcitB, SpinOrbA, OrbA, SymA, SymB
        integer :: nExcitA, sumMl, mlA, mlB, iSpn, Elec1Ind, Elec2Ind
        integer :: SymProduct, ForbiddenOrbs
        logical :: tAOrbFail

!First, we need to pick an unbiased distinct electron pair.
!These have symmetry product SymProduct, and spin pair iSpn = 1=beta/beta; 2=alpha/beta; 3=alpha/alpha
        CALL PickElecPair(nI, Elec1Ind, Elec2Ind, SymProduct, iSpn, SumMl, -1)


!This routine finds the number of orbitals which are allowed by spin,
! but not part of any spatial symmetry allowed unoccupied pairs.
!This number is needed for the correct normalisation of the probability of drawing
! any given A orbital since these can be chucked and redrawn.
        IF (tNoSymGenRandExcits) THEN
            CALL FindNumForbiddenOrbsNoSym(ForbiddenOrbs, ClassCountUnocc2, iSpn)
        ELSE
            CALL FindNumForbiddenOrbs(ForbiddenOrbs, ClassCountUnocc2, SymProduct, iSpn, SumMl)
        end if

!Now we have to pick the first unoccupied orbital. If an orbital is not present in any allowed pairs,
! it is chucked and a new one drawn.
!The number NExcit is the number of unoccupied orbitals that the orbital was chosen from
! (including symmetry-forbidden orbital pairs) Arguments:
        CALL PickAOrb(nI,iSpn,ILUT,ClassCountUnocc2,NExcitA,Elec1Ind,Elec2Ind,SpinOrbA,OrbA,SymA,SymB, &
      &             SymProduct,SumMl,MlA,MlB,ForbiddenOrbs,tAOrbFail)
        IF(tAOrbFail) THEN
            nJ(1)=0
            pGen=HUGE(0.0_dp)
            RETURN
        end if

!This routine will pick an unoccupied orbital at random from a specified spin and symmetry class.
!There should definitely be a possible spinorbital, since A was chosen so that there was one.
!We have to make sure with alpha/alpha or beta/beta pairs and when SymProduct=0,
! that we don't choose the same unoccupied orbital.
!If we do this, then we should chuck and redraw,
! since there should definitely be another allowed spinorbital in the class.
!We return the number of allowed B's for the A we picked in NExcitB,
! however we also need to know the number of allowed A's if we
!had picked B first. This will be returned in NExcitOtherWay.
        CALL PickBOrb(nI, iSpn, ILUT, ClassCountUnocc2, SpinOrbA, OrbA, SymA, OrbB, SymB, NExcitB, MlA, MlB, NExcitOtherWay)

        call make_double(nI, nJ, elec1ind, elec2ind, orbA, orbB, &
                         ExcitMat, tParity)

        CALL FindDoubleProb(ForbiddenOrbs, NExcitA, NExcitB, NExcitOtherWay, pGen)

    END SUBROUTINE CreateDoubExcit