SYMPROD Function

    elemental FUNCTION SYMPROD(ISYM1, ISYM2)
        TYPE(Symmetry), intent(in) :: ISYM1, ISYM2
        TYPE(Symmetry) SYMPROD
        TYPE(Symmetry) IS1, IS2
        INTEGER I, J, Abel1(3), Abel2(3)
        character(*), parameter :: this_routine = 'SYMPROD'
        if (TAbelian) then

            IF (TwoCycleSymGens) THEN
!For molecular systems, we can only have a maximum of 8 irreps, and so can do a simple
!xor to get the symmetry.

                SymProd%s = IEOR(ISym1%s, ISym2%s)

            ELSE

                call DecomposeAbelianSym(ISym1%s, Abel1)
                call DecomposeAbelianSym(ISym2%s, Abel2)
!Slightly faster when calling a lot to do it in an array operation
                Abel1 = modulo(Abel1 + Abel2, NProp)
                SymProd%s = ComposeAbelianSym(Abel1)
            end if
        else
            IF (ISYM1%s == 0 .OR. ISYM2%s == 0) THEN
                SYMPROD%s = 0
                RETURN
            end if
            ! fix if all symmetries are set to 1
!              if (isym1%s == isym2%s .and. isym1%s == 1) then
!                  symprod%s = 1
!                  return
!              end if
            ! change for the real-space hubbard
            IF (.not. allocated(SYMTABLE)) call stop_all(this_routine, 'SYMMETRY TABLE NOT ALLOCATED')
            IS1 = ISYM1
            I = 1
            SYMPROD%s = 0
            if (t_k_space_hubbard) then
                symprod = SYMTABLE(isym1%s, isym2%s)
            else
                DO WHILE (IS1%s /= 0)
                    IF (BTEST(IS1%s, 0)) THEN
                        IS2 = ISYM2
                        J = 1
                        DO WHILE (IS2%s /= 0)
                            IF (BTEST(IS2%s, 0)) THEN
                                SYMPROD%s = IOR(SYMPROD%s, SYMTABLE(I, J)%s)
                            end if
                            !  RSHIFT(,1)
                            IS2%s = ISHFT(IS2%s, -1)
                            J = J + 1
                        end do
                    end if
                    !RSHIFT(,1)
                    IS1%s = ISHFT(IS1%s, -1)
                    I = I + 1
                end do
            end if
        end if
        RETURN
    END FUNCTION SYMPROD