IsMomentumAllowed Function

    FUNCTION IsMomentumAllowed(nJ)

        use sym_mod, only: mompbcsym
        use SystemData, only: kvec, tUEG2, tAllSymSectors

        LOGICAL :: IsMomentumAllowed ! Returns whether the determinant is momentum allowed for
        ! UEG and Hubbard models
        ! Compares the total k from a determinant nI with kTotal
        INTEGER :: nJ(NEl), kx, ky, kz, ktrial(3), i

        if (tAllSymSectors) then
            IsMomentumAllowed = .true.
            return
        end if

        IsMomentumAllowed = .false.

        !====================================================
        if (tUEG2) then
            ! The momentum constraint for UEG: every determinant must have a total momentum
            ! which is equal.
            kx = 0
            ky = 0
            kz = 0
            do i = 1, NEl
                kx = kx + kvec(nJ(i), 1)
                ky = ky + kvec(nJ(i), 2)
                kz = kz + kvec(nJ(i), 3)
            end do
            IF (kx == kTotal(1) .and. ky == kTotal(2) .and. kz == kTotal(3)) THEN
                IsMomentumAllowed = .true.
            end if
            return
        end if
        !====================================================
        ! The momentum constraint for UEG: every determinant must have a total momentum
        ! which is equal.
        IF (tUEG) THEN
            kx = 0
            ky = 0
            kz = 0
            do i = 1, NEl
                kx = kx + G1(nJ(i))%k(1)
                ky = ky + G1(nJ(i))%k(2)
                kz = kz + G1(nJ(i))%k(3)
            end do
            IF (kx == kTotal(1) .and. ky == kTotal(2) .and. kz == kTotal(3)) THEN
                IsMomentumAllowed = .true.
            end if
        end if

        ! The momentum constraint from Hubbard model: every determinant must have a total momentum
        ! which is equal to within a reciprocal lattice vector.
        IF (tHub) THEN
            kx = 0
            ky = 0
            kz = 0
            do i = 1, NEl
                kx = kx + G1(nJ(i))%k(1)
                ky = ky + G1(nJ(i))%k(2)
                kz = kz + G1(nJ(i))%k(3)
            end do
            ktrial = (/kx, ky, 0/)
            CALL MomPbcSym(ktrial, nBasisMax) ! This re-maps the total momentum under PBCs: equivalent to this being equal to
            ! a value to within a reciproval lattice vector.
            IF (ktrial(1) == kTotal(1) .and. ktrial(2) == kTotal(2)) THEN
                IsMomentumAllowed = .true.
            else
                print *, "ERRONEOUS MOMENTUM: ", ktrial, ktotal
            end if
        end if

    END FUNCTION IsMomentumAllowed