GETSYM Subroutine

    SUBROUTINE GETSYM(NI2, NEL, G1, NBASISMAX, ISYM)
        IMPLICIT NONE
        INTEGER NEL, NI(NEL), nBasisMax(5, *)
        TYPE(BasisFn) G1(*), ISym
        INTEGER I, J, NI2(NEL)
        INTEGER NREPS(NEL), NELECS(NEL), SSYM
        I = 1
        NREPS(1:NEL) = 0
        CALL SETUPSYM(ISYM)
        IF (tFixLz) THEN
            CALL GetLz(NI2, NEL, ISYM%Ml)
        ELSE
            ISYM%Ml = 0
        end if
        NI(1:NEL) = NI2(1:NEL)
        IF (tAbelian) THEN !For Abelian symmetry we don't need the symreps malarky.
            DO I = 1, NEL
                ISYM%Sym = SYMPROD(ISYM%Sym, G1(NI(I))%Sym)
!   add the momentum
                CALL ADDELECSYM(NI(I), G1, NBASISMAX, ISYM)
            end do
        ELSE
            DO I = 1, NEL
                !   Count all electrons in each rep
                !   NREPS(J) is the rep, and NELECS(J) is the number of electrons in that rep

                J = 1
                DO WHILE (J < NEL)
                    IF (NREPS(J) == 0) exit
                    IF (NREPS(J) == SYMREPS(1, NI(I))) THEN
                        !   We've found the slot for the rep.  increment it and leave.
                        NELECS(J) = NELECS(J) + 1
                        J = NEL
                    end if
                    J = J + 1
                end do
                IF (J <= NEL) THEN
                    !   need to put the new rep in a new space
                    NREPS(J) = SYMREPS(1, NI(I))
                    NELECS(J) = 1
                end if
            END DO
            !   now go through and see which are closed and which open
            DO I = 1, NEL
                J = 1
                DO WHILE (NREPS(J) /= SYMREPS(1, NI(I)))
                    J = J + 1
                end do
                !   electron NI(I) is in rep NREPS(J)
                IF (NELECS(J) /= SYMREPS(2, NREPS(J))) THEN
                    !   we don't have a closed shell
                    !   add the sym product
!                     print *, "i, ni(i): ", i, ni(i)
!                     print *, "s1,s2: ", isym%sym%s, G1(ni(i))%sym%s
                    ISYM%Sym = SYMPROD(ISYM%Sym, G1(NI(I))%Sym)
!                      print *, "isym after: ", isym%sym%s
                end if
                !   add the momentum
                CALL ADDELECSYM(NI(I), G1, NBASISMAX, ISYM)
            end do
        end if
!   round the momentum
        CALL ROUNDSYM(ISYM, NBASISMAX)
        RETURN
    END SUBROUTINE GETSYM